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Abstract 

We present some recent investigations resulting from the modelling of neural networks as dynamical 
systems, and dealing with the following questions, adressed in the context of specific models. 

(i) . Characterizing the collective dynamics; 

(ii) . Statistical analysis of spikes trains; 

(iii) . Interplay between dynamics and network structure; 

(iv) . Effects of synaptic plasticity. 

*Laboratoire J. A. Dieudonne, U.M.R. C.N.R.S. N 6621, Universite de Nice Sophia-Antipolis 
tlNRIA, 2004 Route des Lucioles, 06902 Sophia-Antipolis, France. 
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The study of neural networks is certainly a prominent example of interdisciplinary research field. 
From biologists, neurophysiologists, pharmacologists, to mathematicians, theoretical physicists, including 
engineers, computer scientists, robot designers, a lot of people with distinct motivations and questions are 
interacting. With maybe a common "Graal" : to understand one day how brain is working. At the present 
stage, and though significant progresses are made regularly, this promised day is however still in a far 
future. But, beyond the comprehension of brain or even of simpler neural systems in less evolved animals, 
there is also the desire to exhibit general mechanisms or principles that could be applied to such artificial 
systems as computers, robots, or "cyborgs" (we think of the promising research field of brain-control of 
artificial prostheses, see for example the web page http://www-sop.inria.fr/demar/indcx_fr.shtml). Again, 
there are many way of tracking these principles or mechanisms. 

One possible strategy is to propose mathematical models of neural activity, at different space and time 
scales, depending on the type of phenomenon under consideration. However, beyond the mere proposal of 
new models, which can rapidly results in a plethora, there is also a need to understand some fundamental 
keys ruling the behaviour of neural networks, and, from this, to extract new ideas that can be tested 
in real experiments. Therefore, there is a need to make a thorough analysis of these models. This can 
be done by numerical investigations, with, very often, the need of inventing clever algorithms to fight 
the hard problem of simulating, in a reasonable time, and with a reasonable accuracy, the tremendous 
number of degree of freedom and the even larger number of parameters that neural networks have. A 
complementary issue relies in developing a mathematical analysis, whenever possible. 

In this spirit, we present in this paper some recent investigations from the authors and his collabora- 
tors, resulting from the modelling of neural networks as dynamical systems. We warn the reader that this 
paper does not intend to be exhaustive and we shall only briefly mention many works which certainly 
would have deserved a longer presentation in a more extensive review: the works by Ermentrout and 
Kopell on phase response theory [67], van Vreeswijk, Sompolinsky and collaborators |177[I178[I176|I175] . 
[120] , Brunei [3__, E_2l EU [72l [___J 1144] , and many others on neural activity, theory of synchronization and 
spike patterns by Seung [181], Bressloff and Coombes [2HI CSIl _29] Timme [H__ □__] Q__6] EE], Jin [87], 
Diesmann [53J are only a few examples of these omissions. 

Beyond the presentation of those results there is also the willing of raising interesting questions 
emerging from this point of view. After a short presentation of neural networks, and how they can be 
indeed modeled as dynamical systems (section[T]), we list 4 of these questions, and address them in specific 
models. 

• Characterizing the collective dynamics of neural networks models. When considering 
neural networks as dynamical systems, a first, natural issue is to ask about the (generic) dynamics 
exhibited by the system when control parameters vary. This is discussed in section [2] 

• Statistical analysis of spikes trains. Neurons respond to excitations or stimuli by finite se- 
quences of spikes (spike trains). Characterizing spike trains statistics is a crucial issue in neuro- 
science. We approach this question considering simple models. This is discussed in section [3J 

• Interplay between dynamics and synaptic network structure. Neural network are highly 
dynamical object and their behavior is the result of a complex interplay between the neurons 
dynamics and the synaptic network structure. In this context, we discuss how the mere analysis 
of synaptic network structure may not be sufficient to analyse such effects as the propagation of a 
signal inside the network. We also present new tools based on linear response theory [151] . useful 
to analysing this interwoven evolution. This is discussed in section 2] 
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• Effects of synaptic plasticity. Synapses evolve according to neurons activity. Addressing the 
effect of synaptic plasticity in neural networks where dynamics is emerging from collective effects 
and where spikes statistics are constrained by this dynamics seems to be of central importance. We 
present recent results in this context. This is discussed in section [5] 

Obviously, the scope of this paper is not to address these questions in a general context. Instead, we 
choose to present simple examples, that one may consider as rather "academic" , for which one can go 
relatively deep, with the idea that such investigations may reveal useful, when transposed to "realistic" 
neural networks. 

1 Neural Networks as dynamical systems 
1.1 From biological neurons and synapses . . . 

A neuron is an excitable cell. Its activity is manifested by local variations (in space and time) of its 
membrane potential, called "action potentials" or "spikes". These variations are due to an exchange of 
ions species (basically Na + ,K + ,Cl~) which move, through the membrane, from the region of highest 
concentration (outside for Na + , inside for K + ) to the region of lowest concentration. This motion does 
not occur spontaneously. It requires the opening/closing of specific gates in specific ionic channels. The 
probability that a gate is open depends on the local membrane potential, whose variations can be elicited 
by local excitations, induced by external currents, or coming from neighbours pieces of membrane (spike 
propagation). Neurons have a spatial structure, depicted in fig [TJ and spikes propagates along this 
structure, from dendrites to soma, and from soma to synapses, along the axon. 




Figure 1: Sketch of the neuron structure. 
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The response of a given neuron to excitations has a wide variability. This variability is not manifested 
by the shape of the action potential, which is relatively constant for a given neuron. Instead, it is revealed 
by the various sequences of spikes a neuron is able to emit. Depending on the excitation, the response 
can be an isolated spike, a periodic spike train, a burst, etc... About twenty different spike trains forms 
are classified in the literature [55] , 

Neurons are connected together. When a spike train is emitted from the soma toward the synapses, 
via the axon, it eventually reaches the synaptic vesicles. Here, a local variation of the membrane potential 
triggers the release of a neurotransmitter into the synaptic cleft. This neurotransmitter reaches by diffu- 
sion the post-synaptic receptors, located on the dendritic spines, and generates a post- synaptic potential 
(PSP). Contrarily to spikes, PSP have an amplitude which depends on the amplitude of the excitation 
and on the synaptic efficacy. Efficacy evolves according to various mechanisms depending on the activity 
of pre- and post-synaptic neurons. Depending on the pre-synaptic neuron and the neurotransmitter used 
by this neuron, the PSP can be either positive or negative. In the first case the pre-synaptic neuron 
and its synaptic connections are called excitatory. Spikes coming from pre-synaptic neuron increase the 
membrane potential of the post-synaptic neuron which is more keen on generating spike trains. Or PSP 
can be negative, corresponding to an inhibitory pre-synaptic neuron. 

Typically, a neuron is connected to many pre-synaptic neurons and receives therefore many excitatory 
or inhibitory signals. The cumulative effects of these signals eventually generate a response of this neu- 
ron's soma that propagates along the axon up to the synaptic tree, then acting on other neurons, and so on. 

From this short description, we can make the following summary. 

• Neurons are connected to each others in a synaptic network with causal (action /reaction) interac- 
tions. 

• Signals exchanged by neurons are spike trains. Spike trains coming from pre-synaptic neurons 
generate a spike train response of the post-synaptic neuron which propagates to other neurons. 

• Spike trains have a wide variability which generates an overwhelming repertoire of collective dy- 
namical responses. 

As an additional level of complexity the structure of the network constituted by synaptic connections 
can also have a wide range of formal, with multiple layers, different species of neurons, etc. Also, a 
very salient property is the capacity that synapses have to evolve and adap^, according to plasticity 
mechanisms. Synaptic plasticity occurs at many levels of organisation and time scales in the nervous 
system [20j . It is of course involved in memory and learning mechanisms, but it also alters excitability 
of brain area and regulates behavioural states (e.g. transition between sleep and wakeful activity). 
On experimental grounds, synaptic changes can be induced by specific simulations conditions defined 
through the firing frequency of pre- and post-synaptic neurons [231 164j , the membrane potential of the 
post-synaptic neuron [10], spike timing [1191 11241 [T9] (see [123] for a review). Different mechanisms have 
been exhibited from the Hebbian's ones [50] to Long Term Potentiation (LTP) and Long Term Depression 
(LTD), and more recently to Spike Time Dependent Plasticity (STDP) [HI HFJ (see [501 [77] [52] for a 
review). 

1 In mathematical models there is no a priori constraint, while in the real world the network structure is constrained by 
genetics. 

2 Note that not only synapses, but also neurons have the capacity of adaptation (intrinsic plasticity 11221 ). We shall not 
discuss this aspect in the present paper. 
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1.2 ... to models 



Regarding the overwhelming richness of behaviors that neuronal networks are able to display, the the- 
oretical (mathematical or numerical) analysis of these systems is at a rather early stage. Nevertheless, 
some significant breakthrough have been made within the last 50 years, as we shall see in a few examples. 
For this, a preliminary modeling/simplification strategy is necessary, that we summarize as follows. 



1.2.1 Fix a model of neuron 



This essentially means: fix an equation or a set of equations describing the evolution of neuron's membrane 
potential, plus, possibly, additional variables (such as the probability of opening/closing ionic gates 
in Hodgkin-Huxley's like models [94]). This choice can be guided by different and, often, mutually 
incompatible constraints. 

• Biological plausibility. 

• Mathematical tractability. 

• Numerical efficiency. 

Regarding the first aspect one may also only focus on a few biological features. Do we want a model 
that reproduces accurately spike shape, or do we simply want to reproduce the variability in spike trains 
responses whereas spike shape is neglected (e.g. represented by a "Dirac" peak) ? Do we want to focus 
on one specific characteristic of spike trains (probability that a neuron fires, probability that two neurons 
fire within a certain time delay....) ? Clearly, there is a large number of neuron models and, as usual, 
models depend on the questions that you ask. Here are a few examples. 



Hodgkin-Huxley model. This model, dating back to 1952 [94], is still one of the best description 
of neuron spike generation and propagation. Thus, it is very good from the point of view of biological 
plausibility. Unfortunately, its mathematical analysis has not been completed yet and it is computational 
time consuming. In this model, the dynamics of a piece of membrane with capacity C m and potential V 
is given by: 

C "^r = -9»»" ,3 ' 1 (^ - EatJ- g K n 4 (V - E K )-g L (V - E L ) + /„ It (1) 
W)t - ^)(l-n)-«K)„==3^p (2, 

1 dm , T , W , x „ / t ,n m°°(V) — ra .. 

W)1I - a m (V)(l-m)-0 m (V)m= ^ (3) 

1 ^ - a h (V ){ l- h )- M V)*= h ^^ (4) 



j(T)dt ' " lx ' r h (V) 

where m, h, n are additional variables, describing the ionic channels activity (see [551 [751 l9"2l 11111 11161 
I134| ). -E/v a , Ex, El are respectively the Nernst potentials of Na + ,K + ions and additional ionic species 
(like Cl~) grouped together in a leak potential E^. gN a ,9K, 9l are the corresponding conductances. 
j(T) is a temperature dependent time scale (equal to 1 at 6.3°C). a, (3 are transitions rates in the 
masters equations (|2l3l4p used to model the transition open/close of ionic channels. Though the complete 
mathematical analysis of this model has not been performed yet, important results can be found in 

[muniois]. 
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Fitzhugh-Nagumo model One can reduce the Hodgkin-Huxley equations in order to obtain an ana- 
lytically tractable model. In this spirit Fitzhugh [70] and independently Nagumo, Arimoto & Yoshizawa 
[133] , considered reductions of the Hodgkin-Huxley model and introduced an analytically tractable two 
variables model 

ef = fx(v,w), 
% = 9x(v,w), 

where e is a small parameter. The index A refers to the control parameters of the system. In the 
FitzHugh-Nagumo model f\(v,w) — v — v 3 — w + I is a cubic polynomial in v and is linear in w, while 
Qx(v, w) — (v — a — bw). The parameters A = (a, b, I) are deduced from the physiological characteristics 
of the neuron. 



Integrate and Fire models Here, one fixes a real number 6, called the firing threshold of the neuron, 
such that if Vfc(t) > 9 then neuron membrane potential is reset instantaneously to some constant reset 
value V rese t and a spike is emitted toward post-synaptic neurons. Below the threshold, < 9, neuron 
fc's dynamics is driven by an equation of form: 

Gk ~dT + gkVk = ik ' ^ 

where CV- is the membrane capacity of neuron k, gk its conductance and ik a current, including various 
term, depending on modeling choices (external current, ionic current, adaptation current). 
In its simplest form equation ([5]) reads: 

dVk _ _Vk + ik_ ^ 
dt Tk Ck 

where gk is a constant, and Tk = ^r- is the characteristic time for membrane potential decay, when no 
current is present. This model has been introduced in |118j . More generally, conductances and currents 
depend on the previous firing times of the pre-synaptic neurons [149] (see section 12.21 for an example) . 



Discrete time models In many papers, researchers use sooner or later numerical simulations to guess 
or validate original results. Most often this corresponds to a time discretisation with standard schemes 
like Euler, or Runge-Kutta. Even when seeking more elaborated schemes such as event based integrations 
schemes [311 1146] . which in principle allows one to handle continuous time, there is in fact a minimal time 
scale, due to numerical round-off error, below which the numerical scheme is not usable anymore. On 
more fundamental grounds, in all models presented above including Hodgkin-Huxley, there is a minimal 
time scale imposed by Physics. Thus, although the mathematical definition of assumes a limit dt — > 0, 
there is a time scale below which the ordinary differential equations lose their meaning. Actually, the 
mere notion of "membrane potential" already assumes an average over microscopic time and space scales. 
Another reason justifying time discretisation in models is the use of "raster plot" to characterize neurons 
activity. 



Raster plots A raster plot is a graph where the activity of a neuron is represented by a mere vertical 
bar each "time" this neuron emits a spike. When focusing on spiking neurons models, spikes are often 
characterized by their "time" of occurrence. Except for IF models, where the notion of "instantaneous" 
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firing and reset leads to nice pathologie^l, a spike has some duration and spike time has some uncertainty 
5. Therefore, the statement "neuron i fires at time t" must be understood as "neuron i fires at time t 
within a precision 5 > 0" . Moreover, a neuron cannot fire more than once within a time period r called 
"refractory period" . Therefore, one can fix a positive time scale S > which can be mathematically 
arbitrary small, such that (i) a neuron can fire at most once between [t, t + #[ (i.e. S << r, the refractory 
period); (ii) dt « 8, so that we can keep the continuous time evolution of membrane potentials, taking 
into account time scales smaller than 5, and integrating membrane potential dynamics on the intervals 
[t, t + S[; (hi) the spike time is known within a precision S (see |f f 4] for an interesting discussion on this 
approach). 

At this stage let us introduce a concept /notation used throughout this paper. One can associate to 
each neuron k a variable u>k(t) = f if neuron k fires between [t, i + and Ufc(i) = otherwise. A "spiking 

pattern" is a vector ui{t) A = [^k{t)]^ =1 which tells us which neurons are firing at time t. A "raster plot" 

is a sequence to = {w(i)}^o, of spiking patterns. We denote [to] Qt — {^(s)}g =0 ' ^ ne ras te r plot from 
time to time t. 

1.2.2 Fix a model of synapse 

Voltage- and activity-based models A single action potential from a pre-synaptic neuron j is seen 
as a post-synaptic potential by a post-synaptic neuron i (see Fig. [I}. The conductance time-course 
after the arrival of a post-synaptic potential is typically given by a function a.y (t — s) where s is the 
time of the spike hitting the synapse and t the time after the spike. (We neglect here the delays due 
to the distance travelled down the axon by the spikes). In voltage-based models one assumes that the 
post-synaptic potential has the same shape no matter which pre-synaptic population caused it, the sign 
and amplitude may vary though [66] . This leads to the relation: 

otij(t) = WijOti(t), 

where cti represents the unweighted shape (called a a-shape) of the post-synaptic potentials. Known 
examples of a-shapes are cti (t) = Kie-^ Ti H(t) or a. t (t) = Kite-^ Ti H(t) where H is the Heaviside 
function. More generally this is a polynomial in t and this is the Green function of a linear differential 
equation of order k: 

±a?^(t)=5(t). (7) 

Wij is the strength of the post-synaptic potentials elicited by neuron j on neuron i (synaptic efficacy or 
"synaptic weight" ) . 

In activity-based models the shape of a PSP depends only on the nature of the pre-synaptic cell, that 
is [66 : 

ctij(t) = Wijdj(t). 

Assuming that the post-synaptic potentials sum linearly, the average membrane potential of neurorQ 
i is: 

3 Consider a loop with two neurons, one excitatory and the other inhibitory. Depending on the synaptic weights, it 
is possible to have the following situation. The first neuron fires instantaneously, excites instantaneously the second one, 
which fires instantaneously and inhibits instantaneously the first, which does not fire... This type of causal paradoxes, 
common in science-fiction novels |15| . can also be found in IF models (eq. 10) without refractory period and time delays. 

4 One should instead write neuron i's soma. In the sequel we shall consider neurons as punctual, without spatial structure. 
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v^) = £M*-4 n) )> (8) 

(n) 

where the sum is taken over the arrival times tj < t of the spikes produced by the neurons j. 

Synaptic plasticity. Most often, the mechanisms involved in synaptic plasticity have been revealed 
by simulation performed in isolated neurons in in vitro conditions. Extrapolating the action of these 
mechanisms to in vivo neural networks requires both a bottom-up and top-down approach. This issue is 
tackled, on theoretical grounds, by inferring "synaptic updates rules" or "learning rules" from biological 
observations [1791 |2"UI I128j (see [501 [TTJ 1^2 for a review) and extrapolating, by theoretical or numerical 
investigations, what are the effects of such synaptic rule on such neural network model. This approach 
relies on the belief that there are "canonical neural models" and "canonical plasticity rules" capturing the 
most essential features of biology. When considering synaptic adaptation, one proposes evolution rules 
for the ctij profiles. Most often, the mere evolution of the Wij's are considered. Here are a few typical 
examples. 

Generic synaptic update Synaptic plasticity corresponds to the evolution of synaptic efficacy which 
evolve in time according to the spikes emitted by the pre- and post- synaptic neuron. In other words, the 
variation of Wij at time t is a function of the spiking sequences of neurons i and j from time t — T s to 
time t, where T s is time scale characterizing the width of the spike trains influencing the synaptic change. 
In its more general form synapse update reads: 

SWij = g (w«(t), {^ } } t , {tf}J , t > T s , 

where {4^} > ({4"^} ) are tnc usts °^ spikes times emitted by the pre-synaptic neuron i, (the post- 
synaptic neuron j), up to time t. Thus, synaptic adaptation results from an integration of spikes over 
the time scale T s . 

With the concept of "raster plot" introduced at the end of section [1.2.21 we may also write: 

SWij = g (Wijit), N*-T s ,t , [u>j]t-T.,t) ^>T S . (9) 
Let us now give a few examples of synaptic adaptation "rules" . 

Hebbian learning In this case, synapses changes depend on the firing rate of neuron i,j. A typical 
example corresponds to 

1 - 

9ij(Wij, N t _ Ts(t , [oJj] t _ T3 it ) = y M s i) - r i(si))(Uj(s 2 ) - r j (s 2 ), (10) 

s Sl:S2= t-T s 

(correlation rule [143] ) where Ti(t) = 4- Y^ s =t~T w i( s ) i s the frequency rate of neuron i in the raster plot 
to, computed in the time window [t — T a ,t], 

5 For further explanations of this terminology, see section [5.21 
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Spike-Time Dependent Plasticity as derived from Bi and Poo [TH] provides the average amount 
of synaptic variation given the delay between the pre- and post-synaptic spike. Thus, "classical" STDP 
reads (TTlfTDU]: 

t 

g(Wij,[Ui]t-T s ,t>{Uj]t-T s ,t) = f(8l~ S2)u}i(si)u}j(s2), (11) 

Sl,S 2 =t-T s 

with: 

{A_e~, x < 0; 
A+e~^, x>0; (12) 
0, x = 0; 

where A- < and A + > 0. The shape of / has been obtained from statistical extrapolations of 
experimental data. Hence STDP is based on a second order statistics (spikes correlations). There is, in 
this case, an evident time scale T s = max(r_, r + ), beyond which / is essentially zero. 
Many other examples can be found in the literature |100] . 

1.2.3 Fix a synaptic graph structure 

This point is closely related to the previous one. In particular, this structure can be fixed or evolve in 
time (synaptic plasticity). In this last case, there is a complex interaction between neuron dynamics 
and synapses dynamics. This structure can be guided from biological/ anatomical data, or it can be 
random. In this last case, one is more interested in generic mathematical properties than by biological 
considerations. The intermediate case can also be considered as well: deterministic synaptic architecture 
with random fluctuations of the synaptic efficacy (see section [2TTI for an example). 

At this stage an interesting issue is : "what is the effect of the synaptic graph structure on neurons 
dynamics ?" This question is closely related to the actual research trend studying dynamical systems 
interacting on complex networks where most studies have focused on the influence of a network structure 
on the global dynamics (for a review, see [H]). In particular, much effort has been devoted to the 
relationships between node synchronization and the classical statistical quantifiers of complex networks 
(degree distribution, average clustering index, mean shortest path, motifs, modularity...) [83 | 1137] 1117] . 
The core idea, that the impact of network topology on global dynamics might be prominent, so that these 
structural statistics may be good indicators of global dynamics, proved however incorrect and some of 
the related studies yielded contradictory results [1371 [53] . Actually, synchronization properties cannot be 
systematically deduced from topology statistics but may be inferred from the spectrum of the network [12] . 
Moreover, most of these studies have considered diffusive coupling between the nodes [89]. In this case, 
the adjacency matrix has real non-negative eigenvalues, and global properties, such as stability of the 
synchronized states [13j can easily be inferred from its spectral properties. 

Unfortunately, this wisdom cannot be easily transposed to the field of neural networks where coupling 
between neurons (synaptic weights) in neural networks is not diffusive, the corresponding matrix is not 
symmetric and may contain positive and negative elements. More generally, as exemplified in sections 
IH and [5] neural networks constitute nice examples where the analysis of the synaptic graph with tools 
coming from the field of "complex networks" provides poor information on dynamics. The main reason 
of this failure is that the synaptic graph does not take into account nonlinear dynamics. In section 2] 
we introduce a different concept of network, based on linear response theory, which provides much more 
information on the conjugated effects of topology and dynamics. 
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1.2.4 Neural networks as dynamical systems 

To summarize, we shall adopt in this paper, the following point of view. "A neural network is for- 
mally a graph where the nodes are the neurons and the edges the synapses, each edge being weighted 
by the corresponding synaptic efficacy. Thus synapses constitute a signed and oriented graph. Each 
node is characterized by an evolution equation where the neuron state depends on its neighbours (pre- 
synaptic neurons). Synaptic weights can be fixed or evolve in time (synaptic plasticity) according to the 
state/history of the two nodes it connects (pre- and post-synaptic neuron)." 

As indicated by the title of this paper we adopt here the point of view that neural networks are 
dynamical systems and we analyse them in this spirit. This point of view is not necessarily completely 
appropriate, but it nevertheless allows some significant insights in neuronal dynamics. More precisely, 
we consider the following setting. 

Canonical formulation of neurons dynamics Each neuron i is characterized by its state, Xj,, which 
belongs to some compact set X S IR M . M is the number of variables characterizing the state of one neuron 
(we assume that all neurons are described by the same number of variables). A typical example is M = 1 
and Xi = Vi is the membrane potential of neuron i and T = [V m i n , V ma x]- Other examples are provided 
by conductances based models of Hodgkin-Huxley type ([T]) then Xi — (Vi,rm,ni,hi) where m^n^ are 
respectively the activation variable for Sodium and Potassium channels and hi is the inactivation variable 
for the Sodium channel. 

We consider the evolution of N neurons, given by a deterministic dynamical system of type: 

<iX 

= F-,, (X, t) , continuous time, (13) 

or, 

X(i + 1) = F 7 [X(t), t] , discrete time. (14) 

The variable X = {Xi}f =1 represents the dynamical state of a network of N neurons at time t. We 
use the notation V instead of X when neuron's state is only determined by membrane potential whereas 
we use the general notation X when additional variables are involved. 

Typically X e M = T N where M is the phase space of (fT4"|) . and F 7 (A4) C M. The map F 7 : M — > 
M depends on a set of parameters 7 € R p . The typical case considered here is 7 = (W,!^^) where 
W is the matrix of synaptic weights and \( ext ) is some external current or stimulus. Thus 7 is a point in 
a P = N 2 + N dimensional space of control parameters. 

Correspondence between membrane potential trajectories and raster plots Typically a neu- 
ron i "fires" (emits a spike or action potential), whenever its state Xi belong to some connected region 
V\ of its phase space. Otherwise, it is quiescent (X e Vo =I\V\). For N identical neurons this leads 
to a "natural partition" V of the product phase space Ai. Call A = {0,1}^, lj — [wjL =1 € A. Then, 
V = {"Pw} w£ a, where — V LUl x x • • • x V UN . Equivalently, if X e V^, then all neurons such that 
u>i = 1 arc firing while neurons such that LOk = are quiescent. 

To each initial condition X 6 M. we associate a "raster plot" u> — such that X(£) € 

^u>(t)>Vi > 0. We write X^u. Thus, ui is the sequence of spiking patterns displayed by the neural 
network when prepared with the initial condition X. On the other way round, we say that an infinite 
sequence u> = {w(i)} t 1 ^ is an admissible raster plot if there exists X G M. such that X — 1 u>. We call E 7 
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the set of admissible raster plots for the set of parameters 7. The dynamics of X induces a dynamics on 
the set of raster plot given by the left shit er 7 such that a~,uj = u>' u'{t) = u(t + l),Vt > 0. Thus, in 
some sense, raster plots provide a code for the orbits of ([Tl| . Note that the correspondence may not be 
one-to-one. 

2 Collective dynamics 

When considering neural networks as dynamical systems, a first, natural issue is to ask about the (generic) 
dynamics exhibited by the system when control parameters (summarised by the symbol 7 in the section 
ll.2.4p vary. However, at the present stage, this question is essentially unsolvable, taking into account 
the very large number of degree of freedom and the even larger number of parameters. Also, the mere 
notion of genericity has to be clarified. In dynamical systems theory "generic" has two distinct meanings. 
Either one is seeking properties holding in a residuajj set, in which case one deals with genericity in a 
topological sense. Or one is interested in properties holding on a set of parameters having probability 
one, for a smooth and "natural" probability distribution defined on the space of control parameters (e.g. 
Lebesgue or Gauss distribution). In this case, one speaks about "probabilistic genericity". These two 
notions of genericity usually do not coincide. (An attempt to unifying these two concepts has been 
proposed in [98] under the name of "prevalence" ) . 

Genericity results are relatively seldom in the field of neural networks, unless considering some specific 
situations (e.g. weakly coupled neural networks, where some neurons of the uncoupled system, are close 
to the same codimension one bifurcation point [95] )■ We present here two genericity results in this section, 
and the related techniques. For a wider review see [1541141) . See also [166] for a new and recent approach. 

2.1 Mean-field methods. 

As a first example let us describe within details the so-called dynamic mean-field theory. This method, 
well known in the field of statistical physics and quantum field theory, is used in the field of neural 
networks dynamics with the aim of modeling neural activity at scales integrating the effect of thousands 
of neurons. This is of central importance for several reasons. First, most imaging techniques are not 
able to measure individual neuron activity ("microscopic" scale), but are instead measuring mesoscopic 
effects resulting from the activity of several hundreds to several hundreds of thousands of neurons. Second, 
anatomical data recorded in the cortex reveal the existence of structures, such as cortical column^], with 
a diameter of about 50/ito to 1mm, containing of the order of one hundred to one thousand neurons 
belonging to a few different species. In this case, information processing does not occur at the scale 
of individual neurons but rather corresponds to an activity integrating the collective dynamics of many 
interacting neurons and resulting in a mesoscopic signal. 

6 A set is residual if is the countable intersection of open dense sets. In this context, "generic" means "holding on a dense 
set of parameters" . 

7 Cortical columns are small cylinders, of diameter ~ 0.1 — 1 mm, that cross transversely cortex layers. They are involved 
in elementary sensory-motor functions such as vision. They are composed of several hundred to thousand neurons, belonging 
to a few different populations belonging to distinct cortex layers. The electrical activity of cortical columns can be measured 
using different techniques. In Optical Imaging, one uses Voltage-Sensitive Dyes (VSDs). The dye molecules act as molecular 
transducer that transform changes in membrane potential into optical signals with a high temporal resolution, < 1 ms, 
and a high spatial resolution, ~ 50 fim. The measured optical signal is locally proportional to the membrane potential of 
all neuronal components and proportional to the excited membrane surface of all neuronal components I81| . It is possible 
to propose phenomenological models characterising the mesoscopic electrical activity of cortical columns. This is useful to 
predict the behaviour of the local field potential generated by neurons activity and to compare this behaviour to measures. 
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However, obtaining the equations of evolution of the effective mean-field from microscopic dynamics 
is far from being evident. In simple physical models this can be achieved via the law of large numbers and 
the central limit theorem, provided that time correlations decrease sufficiently fast. The idea of applying 
mean- field methods coming from statistical physics to neural networks dates back to Amari pj [5] . Later 
on, Crisanti, Sompolinsky and coworkers [163] used a dynamic mean-field approach to conjecture the 
existence of chaos in an homogeneous neural network with random independent synaptic weights. This 
approach was formerly developed by Sompolinsky and coworkers for spin-glasses |164[ 1571 156j . Later 
on, the mean-field equations derived by Sompolinsky and Zippelius [164] for spin-glasses were rigorously 
obtained by Ben Arous and Guionnet [TH] El [53] • The application of their method to a discrete time 
version of the neural network considered in [163j and in [131] was done by Moynot and Samuelides . 132 , . 
Alternative approaches have been used to get a mean-field description of a given neural network and to 
find its solutions. A static mean-field study of multi-population network activity was developed by Treves 
in [173] . His analysis was completed in [1] , where the authors considered a unique population of nonlinear 
oscillators subject to a noisy input current. They proved, using a stationary Fokker-Planck formalism, 
the stability of an asynchronous state in the network. Later on, Gerstner in [75] built a new approach 
to characterize the mean-field dynamics for the Spike Response Model, via the introduction of suitable 
kernels propagating the collective activity of a neural population in time. Brunei and Hakim considered 
a network composed of integrate-and-fire neurons connected with constant synaptic weights [32]. In the 
case of sparse connectivity, stationarity, and considering a regime where individual neurons emit spikes at 
low rate, they were able to study analytically the dynamics of the network and to show that the network 
exhibits a sharp transition between a stationary regime and a regime of fast collective oscillations weakly 
synchronized. Their approach was based on a perturbative analysis of the Fokker-Planck equation. A 
similar formalism was used in [125j which, when complemented with self-consistency equations, resulted 
in the dynamical description of the mean-field equations of the network, and was extended to a multi- 
population network. Finally, Chizhov and Graham [48j have recently proposed a new method, based on 
a population density approach, allowing to characterize the mesoscopic behaviour of neuron populations 
in conductance-based models. 

The motivations of this section are twofold. On the one hand, we present an example of dynamic mean- 
field approach applied to plausible models of mesoscopic neural structures in the brain [69] . Especially, we 
insist on the rich phenomenology brought by this method. On the other hand we present some examples of 
bifurcations analysis of dynamical mean-field equations and what this tells us about the generic dynamics 
of the underlying neural network. 



2.1.1 Multi-populations dynamics 

Brain structures such as cortical columns are made of several species of neurons (with different physical 
and biological characteristics) linked together in a specific architecture |168j . We model this in the 
following way. We consider a network composed of N neurons indexed by i € {1, . . . , N} belonging to 
P populations indexed by a £ {1, . . . , P}. Let N a be the number of neurons in population a. We have 
N = E«=i Na - In thc following we are interested in the limit N —> oo. We assume that the proportions 
of neurons in each population are non-trivial, i.e. : 

lim ^= Pa €]0,1] ;Vae {1, P}. 

N— >oo iV 

On the opposite, were p a to vanish, would the corresponding population not affect the global behavior of 
the system and could it be neglected. We introduce the function p : {1, . . . , N} — * {1, . . . , P} such that 
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p(i) is the index of the population which the neuron i belongs to. 



Firing rates models. In many examples the spiking activity is resumed by spike rates. Call Vj(t) the 
spikes rate of neuron j at time t such that the number of spikes arriving between t and t + dt is Vj(t)dt. 
Moreover, the relation between the membrane potential of neuron i, Vi and Vi takes the form: 



Vi (t) = Si{Vi{t)), 

[77l [60] , where Si is sigmoidal. Therefore, we have, for voltage-based models, 



V, 



(t) = £ a t {t-s) Wi^iVjis)) + Iiis) + Bi(s) ds 



J2w t3 S J (V 3 ) + I l +B l 



(15) 



(f), (16) 



where * is the convolution product. Here we have assumed that neuron i receives also an external current 
Ii(t) and (white) noise Bi(t). For activity based models, defining the activity as: 



one has 



Aj(t) = / aij(t- s)vj(s)ds, 



+ a % * Ii + a z * Bi 



(17) 



Model dynamics Applying the Green relation ([7]) to the membrane potential of the voltage based 
model ([TBI) one obtains: 



N 



1=0 



dt 1 



The first term of the l.h.s. is the contribution of the pre-synaptic neurons to the time variation of 
the membrane potential. Under the assumption that the a-shape, sigmoidal shape, external current and 
noise only depend only on the neuron's population we may write, for each neuron in the population a: 



,d l Vi 



P N b 



E o - -*r(*) = EE^WKt)) + i a {t) + B a (t), i e a. 



dt 1 

1=0 6=1 3=1 

In the case where a a = e~~ H(t) p8|) becomes: 



(18) 



dV V- P Nh 

= -f + E E WiiSM®) + + B S), i e a. 

Ta 6=1 j=l 



(19) 



called the "simple model" in the sequel. 
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Synaptic weights When investigating the structure of mesoscopic neural assemblies such as cortical 
columns, experimentalists are able to provide the average value of the synaptic efficacy from a neural 
population to another one 168J. Obviously, these values are submitted to some indeterminacy (error 
bars). We model this situation in the following way. Each synaptic weight Wij is modeled as a Gaussian 
random variable whose mean and variance depend only on the population pair a = p(i), b = p(j), and on 
the total number of neurons Nb of population b: 



2 
6 

N b ' N b 



where Af(m,a) denotes the Gaussian law with mean m and variance a. We assume that the Wij's are 
uncorrelated. We use the convention Wij — whenever there is no synaptic connection from j to i. 

2.1.2 Mean-Field approach 

Local interaction field The collective behaviour of neurons in eq. (fT8|) is determined by the term: 

P N b 
6=1 j=l 

called the "local interaction field" of neuron %. When the W^'s are fixed, its evolution depends on the 
evolution of all neurons (i.e. the trajectory of the corresponding dynamical system). If the trajectory is 
prescribed, and if the Wij 's vary, rji (i) becomes a random process whose law is constrained by the law of the 
's. Let us analyse this within more details. We first make a qualitative description explaining the basic 
ideas without mathematical rigor. Especially, we assume that there is a well defined "thermodynamic 
limit" (N — > oo) for the quantities we consider. Then we quote a rigorous result validating this qualitative 
description [69]. It uses large deviations techniques developed in [T6 l fT7 l [86 lll32j (see |154j for a review). 

Non random synaptic weights Assume that a a b — 0, namely we neglect the errors in the synaptic 
weights determination. Then, we may write: 

p W Nb 

%(*) = E^E^(*))- ( 2 °) 

6=1 j=l 

As Ni, — > oo, 

1 N b 

— J2S b (Vj(t))^WV(t)), (21) 
j=i 

assuming that the limit exists. The quantity 0b(V(£)) is the average firing rate of population b at time 
t. In this limit, eq. (|18p becomes: 

k i P 

E (*) = E WMY{t)) + I a (t) + B a (t), i e a. 

1=0 6=1 



In this equation the membrane potential evolution only depends on the neuron's i population. Thus, 
l 



setting V a (t) = lim^—oo Y^iti we nave: 
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E ffl i !) ^rW = E W^faMt)) + '«(*) + B a (t), a=l...P, (22) 

(=0 b=l 

called the "first order mean-field" equations in the sequel. 

This equation resembles very much eq. (|18[) if one makes the following reasoning. "Since 0f,(V(i)) is 
the frequency rate of neurons in population 6, averaged over this population, and since, for one neuron, 
the frequency rate is Vi(t) = Si(Vi(t)), let us write 



&(V(t)) = S b (V b (t)). 

This leads to: 

k 1 P 

E °£ ^r(*) = E W ab S b (V b (t)) + I a (t) + B a (t), a = 1 . . . P, 

(=0 b=l 

which has exactly the same form as eq. (fTE)) but at the level of a neurons population. Equations of 
this type, called "naive mean-field" equations in the sequel, are therefore obtained via a "questionable" 
assumption: 

— J2s b (v J (t)) = s b 

b J=l 

There are many examples in physics where this assumption is wrong (such as spin-glasses). However, in 
the present context where the Wij's are independent (and in particular non symmetric, contrarily to e.g. 
spin glasses |127j ) it is correct in some specific sense, as we develop. Actually, naive mean- field equations 
are commonly used as phenomenological models in the neuroscience literature. Here is an example. 




The Jansen-Rit model cortical columns model |102j features a population of pyramidal neurons 
that receives excitatory and inhibitory feedback from local inter-neurons and an excitatory input from 
neighboring cortical units and sub-cortical structures such as the thalamus (see Fig. [5]). The excitatory 
input is represented by an external stimulus with a deterministic part I(t), accounting for some specific 
activity in other cortical units, and a stochastic part B(t) accounting for non specific background activity. 

Denote by V (resp £, I) the pyramidal (respectively excitatory, inhibitory) populations. Choose 
in population V (respectively populations £, T) a particular pyramidal neuron (respectively excitatory, 
inhibitory inter-neuron) indexed by i pyr (respectively i e xc, iinh)- The equations of their activity variable 
read, in agreement with (fl7|) : 

' A lpyr = a £ * S(E WijAj + E WijAj + a £ * /(•) + a £ * £(•)) 

Jexc Jinh 

A iinh =a x *S{Y J W ij A j ) 

Jpyr 

This is therefore an activity-based model. The transfer functions a £ and ax correspond respectively 
to excitatory and inhibitory post-synaptic potential (EPSP or IPSP). In the model introduced originally 



1G 




Figure 2: Schematic representation of the neural populations and their interactions, as considered in Jansen-Rit's model 



by Jansen and Rit, the synaptic integration is of first-order a(t) — Ke~i H(t), where the coefficient K, r 
are the same for the pyramidal and the excitatory population (denote them by K £ ,t £ ), and different 
from the ones of the inhibitory population (denote them by Kx, ri). The sigmoid functions are the same 
whatever the populations. In Jansen-Rit's approach the connectivity weights are assumed to be constant, 
equal to their mean value. Their equations, based on a naive mean-field approach, read therefore, with 
our notations [1021152] : 

^(*) =-TT + K £ S(W V£ A £ (t) + W vx Aj(t) + a £ * /(•) + a £ * £(•)), 
^(i) =-M +K£ s{WsvMt)), (23) 
^(t) =-% + K x S(W X pA v (t)). 

A higher order model, where a(t) — K^e~i H(t), was introduced by van Rotterdam and colleagues 
[174] to better account for the synaptic integration and to better reproduce the characteristics of real 
EPSP's and IPSP's. The bifurcation diagram of this version is quite richer than the Jansen- Rit one [82] , 
These equations are currently used in the neuroscience community either to provide activity models used 
for the analysis of signals obtained from imaging (MEG or Optical Imaging), or to provide dynamical 
models of epilepsy [53] . 

Role of synaptic weights variability Let us now consider the more general case where synaptic 
weights have fluctuations about the mean value W a b- These variations dynamically differentiate the 
neurons within a population and may induce dramatic collective effects, when amplified by the nonlinear 
dynamics. Then, the actual evolution of a population can depart strongly from the first order mean-field 
approximation (not to speak of the naive mean-field approach) . 

Consider the local interaction field (12"0|) . Fix the trajectory of V = {Vi}- =1 . Then, the Wij's being 
Gaussian, rji(t) is (conditionally) Gaussian, with mean: 

p ] N b 

e[«(*)iv]=i;^ f i;w))- 

6=1 j=l 

where E [] is the expectation with respect to the Wij's distribution , and covariance: 
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P 2 N b 

cov Mtfowv] = % E#E W(t))s 6 (y,(*)), 

6=1 j = l 

where we have used that the Wy's are independent so that Cov(Wij, Wi'f) — o~ijO~i'j>j^, i £ a, j £ b. 
Thus, conditionally to V, and still assuming that there is a well defined thermodynamic limit, r\ a converges 
as N — > oo to a diagonal Gaussian process r\ a whose law depends only on the population, with mean: 



E [r)a{t)\V] = WabMV(t)), (24) 

and covariance: 



6=1 



p 1 N b 

Cov [ V a(t)va(s)\V] = J2 ff a6 lim — SbWiityS^Vjia)) (25) 

Thus for a fixed trajectory, we find that the average value of r\ a obeys the same equation as in the 
first order mean-field approach, but it has now fluctuations and correlations given by 



The main difficulty is obviously that the trajectory V is generated by dynamics including the nonlinear 
and collective effects summarized in r\ a . The following result can be shown [53] ■ As N a — > oo the 
membrane potential of a neuron in population V a obeys the equation: 



k H'V P 

'■a 

1=0 6=1 



ft 



^ a alt 1 



(t)=J2Uab(t)+Ia{t)+B a (t) (26) 



where U ab , called the "mean-field interaction process", is a Gaussian process, (thus entirely defined by 
its mean and covariance), statistically independent of the external noise B and of the initial condition 
V(i ), and defined by: 



E [U ab (t)\ = W ab m b {t) where m b (t) d = E[S b (V b (t))}; 
Cov(Uab(t), Uab(s)) = a 2 ab A ab (t,s) where 
A ab (t,s) d = m[s b (V b (t))S b {V b (s)) 
{Gov{U ab (t), U cd (s)) = if ft ^ c or d. 



(27) 



One obtains therefore a set of self-consistent equations giving the mean and covariance of the mean- 
field interaction process U ab . The interaction field of population a, rj a , is given by rj a — U ab , so that 
•q a is indeed a Gaussian process with mean X)&=i W ab mb(t) in agreement with eq. (|24|) . and covariance 
S&=i a ab^ab{t, s) , in agreement with eq. (|2"5|) . But there is a important distinction. Eq. (|2l)j) . (f2"7) 
provide the law of U ab and 7y a , and provide a closed system of equations ruling the dynamical evolution 
of V a averaged over the distribution of W^'s, while in equations (j2"4"|) , (j2"5j) we only got the conditional law 
with respect to a fixed trajectory V, henceforth leading to an incomplete formulation of the problem, 
since, to close the equations, one needs to know the probability distribution of the trajectories V. This 
is an important distinction explaining the difference of notation between </>b(V(£)) in cq. (|24|) and m b (t) 
in mh. 
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Example: the simple model Since U ab is a Gaussian process it is straightforward to obtain an explicit 
form for its mean and covariance as well as for the mean and covariance of V a . In the case of the simple 
model (eq. dHJ)) this leads to the following equation for the evolution of the average value n a (t) of V a : 



dHa _ Ma 
dt t. 



with: 



- + Y. W " b l + S b (hy/^ + to(tj)Dh + I a (t), (28) 

a J-oo V ' 



Dh = ^jJLdh, (29) 

V 27T 



where v a {t) is the variance of V a at time t. Let C a b(t,s) be the covariance of V a (t),Vb(s). Then, 
v a {t) = C ab (t,t). C a b{t, s) is given by [69]: 



C„ 6 (M) = <Wr (t+s)/T ' 
where: 

A b (u, v) = S b I x 



v a (0) + ^ (e£ - l) + [ S e^+^A b (u,v)dudv 

^ •/ «/ 



(30) 



^/ t>b(u)ffc(^) - C b(M, t;) 2 C bb (u,v ) 
\/v b {v) \/vb{v) 



+ V ^j== + fib(u) ) S b (yy/v b (v) + Hb(v)) 
\/v b (v) / \ 7 



Dx Dy, (31) 



and where s 2 is the variance of a white noise Bi(t) in (|16p and where Dx, Dy are Gaussian integrands of 
type ([29]). 

These equations extend as well to more complex models, including the cortical columns model of 
Jansen-Rit [102] and van Rotterdam and colleagues [174] (see [69]). Therefore, the introduction of 
fluctuations in the synaptic distribution change drastically equations of evolution of such neural masses 
models as Jansen-Rit (see [69j for further comments). 



2.1.3 Bifurcations of mean-field equations: a simple but non trivial example 

Let us investigate these equations within details. In the case where fluctuations are neglected (cr 2 ^ = 
0, s a = 0), equations (l30l) .(|3i p admit the trivial solution C ab (t, s) = 0, v b (t) = and equation (|28f reduces 
to the equation obtained by the naive mean-field approach. Incidentally, this validates the naive mean- 
field approach in this context. However, as soon as a\ h > dynamics become highly non trivial since 
the mean-field evolution (|28|) depends on its fluctuations via the variance v b (t). This variance is in turn 
given by a complex equation requiring an integration on the whole past. Actually, unless one assumes 
the stationarity of the process, this equation cannot be written as an ordinary differential equation and 
the evolution is non-Markovian. This result, well known in the field of spin-glasses [18] , has only been 
revealed recently in the field of neural masses models 69J, though mean- field approaches were formerly 
used |163[ [39] [36] . In these last papers, the role of mean-field fluctuations was clearly revealed and its 
influence on dynamics emphasized. In particular, chaotic dynamics have been exhibited, while the mean 
value Ha(t) has a very regular and non chaotic behaviour (for example, it can be constant). 
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The model As an example, let us consider the following model, corresponding to a time discretisation 
of (|19p with dt = t and only one population. Thus, synaptic weights are Gaussian with mean and 

2 

variance jr. Dynamics is given by: 

N 

V i {t + l)=Y / W ij S(V j (t))+I i , i = l..N, (32) 
j'=i 

where S is a sigmoidal function such as S(x) — tanh(g:z;) or S(x) — 1 + ta ° h (g ;r ) . xhe parameter g controls 
the non-linearity of 5*. There is a time-constant current I whose components U are random variables 
with mean I and variance erf. 

The mean-field equations. They write [39|. 136] : 

H{t + 1) = W [ S(h^/v(fj + n(t))Dh + 1, (33) 
Jr 

v(t + l) = o 2 ( S 2 (hy/^t)+ fj,(t))Dh + aj, (34) 
Jr 

C(t+l,t'+l) = a 2 [ s( V<t)v(t') Z C 2 ( tJ) h + 9Ml h i + ^ t )\ S ( h'^Af) + W)) DhDh'+aj, 

Jr 2 V vM^) V^') J v ' 

(35) 

where we have made v(t) explicit, though it can be obtained from (1351) . 

Let us comment these equations. First, they contain statistical parameters determining the probability 
distribution of synaptic weights and currents, W , c, / , err. They also contain an hidden parameter, g 
determining the gain of the sigmoid, which is the same for all neurons. As we saw, deriving mean-field 
equations corresponds to substituting the analysis of the dynamical system (|14[) , with a huge number of 
random microscopic parameters, by an "averaged" dynamical system depending on these few deterministic 
macroscopic parameters. In this spirit, we expect these equations to give indications about the generic 
behavior (in a probabilistic sense) when the synaptic weights and couplings are drawn according these 
values of macroscopic parameters, and when the number of neurons is large. 

The variables m, C essentially play the role of order parameters in statistical physics. They character- 
ize the emergent behavior of a system with a large number of degree of freedom and they exhibit drastic 
changes corresponding, in statistical physics, to phase transitions, and in our context to a macroscopic 
bifurcations. 

Bifurcations in mean-field equations. Having these equations in hand, the idea is now to study the 
reduced dynamical system (J3"3")) . (f3"4"|) , (|35[) and to infer information about the typical dynamics of (f3"2"]) . In 
the present example there exists a stationary regime of ([53")) . (|3"4"|) . (|55|) and the stationary equations are 
given by: 

H = W I S(hy/v + n)Dh + I, (36) 
Jr 

v = a 2 [ S 2 (hy^ + ^)Dh + aj, (37) 
Jr 
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C(t - t') = a 2 [ S 



^yi~C 2 {t-t>) h | Cjt-t') 




fc' + Js (h'y/v + (i) DhDti + aj. 



(38) 



Jr 2 



These equations give important information about the statistical behavior of the model (|3"2"|) with an 
increasing accuracy when the size increases. For example saddle-node bifurcations can be exhibited 
giving rise to bi-stability (see fig. [3] and [39] for more details). But, the most salient feature, as revealed 
by a detailed analysis of the complete set of equations (|36|) . (|37p . ([38| . and especially of the equation 
for the time covariance (|38p). is the existence of a chaotic regime, occurring for a sufficiently large non- 
linearity g. This chaotic region is delimited, in the space of parameters (g, W, a, I, 07) , by a manifold 
whose equation is: 



Note that, in the case where S(x) = ta,nh(gx) this equation gives precisely the so-called De Almeida 
Thouless line [HI], delimiting, in the Sherrington-Kirckpatrick model of spin-glasses [158] . a frontier in 
the plane temperature-local external field, below which dynamics becomes highly non trivial. Here the 
parameter corresponding to the inverse temperature is g, while the external local field corresponds to I. 
The "low temperature regime" of the SK model corresponds therefore to the chaotic regime of (|3"2"|) . This 
analogy is further discussed in J35[ [36] . 

Interpretation. It can be shown that the crossing of this manifold corresponds, in the infinite system, 
to a sharp transition from fixed point to infinite dimensional chaos [35, 36J. Considering the finite size 
system, one can show that ([32]) exhibits generically a Ruelle Takens [153] transition to chaos as g increases. 
As N increases the transition to chaos occurs on a g range becoming more and more narrow, giving this 
sharp transition in the thermodynamic limit. This is related to the fact that the eigenvalues of the 
Jacobian matrix accumulate on the stability circle as N — > oo [78] (see [36j [41] for details) . 

The interesting remark is that, considering only the naive mean-field equation (equation for the 
mean fj,(t) with variance v(t) = 0), one can easily exhibit examples (e.g. S(x) — taiih(gx) with no 
current) where n(t) is a constant (0), while fluctuations are chaotic. This clearly shows the limits of 
the naive mean-field approach and the interest of analysing the role of fluctuations, not only in simple 
models such as ([32]) but also for more realistic models with several populations, like Jansen-Rit's ([23]) . 
Field fluctuations could reveal effects that do not appear in the naive mean-field approach and that 
could be measured in experiments. This question is under investigations (see the web page http:/ /www- 
sop. inria.fr/odyssee/contracts/MACACC/macacc. html for more details). 

2.2 Dynamics of conductance based Integrate and Fire Models 

Let us now investigate a second type of collective dynamics, in the context of Integrate and Fire models 
introduced in section [T. 2. 11 These models have known a great success due to their (apparent) conceptual 
simplicity and analytical tractability 130. 68, 157, 169, 126, IT9"l["101| that can be used to explore some gen- 
eral principles of neurodynamics and coding. Surprisingly, the analysis of only one IF neuron submitted 
to a periodic current reveals already an astonishing complexity and the mathematical analysis requires 
elaborated methods from dynamical systems theory [110, 50, 5Tj. In the same way, the computation of 
the spike train probability distribution resulting from the action of a Brownian noise on an IF neuron 
is not a completely straightforward exercise [1 151 [77l [34l [33l I172j and may require rather elaborated 
mathematics. At the level of networks the situation is even more complex, and the techniques used for 




(39) 
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Figure 3: Schematic bifurcation map of model l|32| l in the plane I,g (drawn from [39]). In the insets are represented the 
type of bifurcation occuring, for the mean-field equations, when crossing the line indicated by an arrow while increasing 
g and keeping I constant. The bifurcation occurs for a value denoted go in the inset. In these insets, the coordinates are 
(g, fi) where fi is given by eq. {36}. The transition from a stable fixed point to chaos corresponds, for finite size systems to 
a Ruelle-Takens transition 153 , while in infinite dimension (TV — » oo) this is a sharp transition (infinitely many eigenvalues 
of the Jacobian matrix at the fixed point crossing simultaneously the frontier of stability). Region I corresponds to a regime 
with one stable fixed point; region II to the coexistence of two stable fixed points; region III to a chaotic attractor and 
region IV to the coexistence of a stable fixed point and of a chaotic attractor. 



the analysis of a single neuron are not easily extensible to the network case. For example, Bressloff and 
Coombes [35J have extended the analysis in [1101 I5U1 I5T] to the dynamics of strongly coupled spiking 
neurons, but restricted to networks with specific architectures and under restrictive assumptions on the 
firing times. Chow and Kopell |49j studied IF neurons coupled with gap junctions but the analysis for 
large networks assumes constant synaptic weights. Brunei and Hakim [33] extended the Fokker-Planck 
analysis combined to a mean-field approach to the case of a network with inhibitory synaptic couplings 
but under the assumptions that all synaptic weights are equal. However, synaptic weights variability 
plays a crucial role in the dynamics, as we saw in the previous section (see also [1761 1178[ 1175] ). Note 
that the rigorous derivation of the mean- field equations, that requires large-deviations techniques [18] . 
has not been yet done for the case of IF networks with continuous time dynamics (for the discrete time 
case see [1651 1154] V 

In this section, we present a rigorous result characterising the generic dynamics of a Generalised 
Integrate and Fire model, where time has been discretized according to the discussion of paragraph 
"raster plots" in section fl. 2. II We then give an example where we consider random synaptic weights. 
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2.2.1 Model 



As we saw, the occurrence of a post-synaptic potential on synapse j, at time vp , results in a change 
of membrane potential (eq. ©). In conductance based models |149] this change is incorporated in the 
adaptation of conductances. The evolution of V k , the membrane potential of neuron k, reads, setting the 
membrane capacity C k = 1 for simplicity: 

^ = (Vk ~ E L ) - i^ n) {V k ,t, M 0)t ) + *<"*>(*), (40) 
at tl 

where [ui] t is the raster plot up to time t. Recall that knowing [oj] t is equivalent to knowing the list 
ji^ j of firing times of all neurons up to time t. The first term in the r.h.s. is a leak term, i k ext \t) is 
an external current, while: 



je£ jei 



" {V k ,t, H 0>t ) = (V k - E+) J29k 3 (t, [w] 0it ) + (V k - E-) J29kj(t, [u\ Q J, 



where E^ are reversal potential (typically E + ~ OmV^ and i£ _ ~ — 75m V). As in the previous section, £ 
andX refers respectively to excitatory and inhibitory neurons, and the + (— ) sign is relative to excitatory 
(inhibitory) synapses. Note that conductances are always positive thus the sign of the post-synaptic 
potential is determined by the reversal potentials E ± . At rest (V k ~ — 70mV) the + term leads to a 
positive PSP while — leads to a negative PSP. 

Conductances depend on past spikes via the relation: 

M 3 -(t,V) 
n=l 

In this equation, Mj(t, V) = ^1=0^ ( s ) ^ s the number of times neuron j has fired at time t. G k j is a 
positive constant proportional to the synaptic efficacy 

W kj = E+G k3 if j G £, 
W fej = £"G fej if j G I. 

Recall that we use the convention W k j = if there is no synapse from j to k 
Then, we may write (|40j) in the form : 

dV k . 
—jj- + g k V k = i k , 

(eq. ^ introduced in section [T. 2. ip with: 

1 W 

fffc(*. Mo,t) = — +^29kj(t, [uj} o t ), 
i=i 

and: 

[«] 0)t ) = ^ + £+ M ,t) + £~ X>i(*> m 0i4 ) + 4 e " t} c*)- 
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This equation characterises the membrane potential evolution below the threshold 6. Recall that, in 
Integrate and Fire models, if V k (t) > then neuron membrane potential is reset instantaneously to some 
constant reset value V rese t and a spike is emitted toward post-synaptic neurons. 

2.2.2 Time discretisation 

Using a time discretisation with a time step 5 = 1, with the hypothesis discussed in section [1.2.11 leads 
to the following discrete-time model |45| : 

V k (t + 1) = 7fc (t, [a;] 0>t ) [1 - Z(V k (t))] V k (t) + J k (t, H 0it ), (41) 

where: 

7*(*,H 0lt ) = e--^ 1 *W w W)«"<i, (42) 
is the integrated conductance over the time interval [t, t + 1[, 

/•t+l 

M ,t) = J ik(s, [w] ,t) ^( s > i + !' Ho,t) ds > 
is the corresponding integrated current with: 

v k (s,t+l,[u] , t ) = e-I« +1 3 ^ s ''M°.J ds ', 

and where Z is defined by : 

Z(x) = X [x > 0} , (43) 

where x lS the indicator function that will later on allows us to include the firing condition in the evolution 
equation of the membrane potential (see (fl4|) ). 

2.2.3 Generic dynamics 

It can be shown that this systems has the following properties. 

Singularity set. The dynamics (|14l) (and the dynamics of continuous time IF models as well) is not 
smooth, but has singularities, due to the sharp threshold definition in neurons firing. The singularity set 
is: 

S = {V G M\Bi = 1 . . . N, such that V t = 8} . 

This is the set of membrane potential vectors such that at least one of the neurons has a membrane 
potential exactly equal to the threshold H. This set has a simple structure: it is a finite union of N — 1 
dimensional hyperplanes. Although S is a "small" set both from the topological (non residual set) and 
probabilistic (zero Lebesgue measure) point of view, it has an important effect on the dynamics. 

8 A sufficient condition for a neuron i to fire at time t is Vi(t) = 8 hence V(t) G S. But this is not a necessary condition. 
Indeed, there may exist discontinuous jumps in the dynamics, even if time is continuous, either due to noise, or a profiles 
with jumps (e.g. a(t) = Ke~~ , t > 0). Thus neuron i can fire with Vi(t) > and V(t) ^ S. In the present case, this 
situation arises because time is discrete and one can have V(t — 5) < 8 and V(t) > 8. This holds as well even if one uses 
numerical schemes using interpolations to locate more precisely the spike time |88| . 
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Local contraction. The other important aspect is that the dynamics is locally contracting, because 
7fc(t, Mot) < 1 ( see ec l- 6210 ■ This nas * ne following consequence. Let us consider the trajectory of a 
point V £ M and perturbations with an amplitude < e about V (this can be some fluctuation in the 
current, or some additional noise, but it can also be some error due to a numerical implementation). 
Equivalently, consider the evolution of the e-ball B(V, e). If B(V, e) n S = then the image of B(V, e) is 
a ball with a smaller diameter. This means, that, under the condition B(V, e) n S = 0, a perturbation is 
damped. Now, if the images of the ball under the dynamics never intersect <S, any e-perturbation around V 
is exponentially damped and the perturbed trajectories about V become asymptotically indistinguishable 
from the trajectory of V. This means that, if the membrane potential of neurons do not approach the 
threshold within a distance smallei@than e then perturbations of size smaller than e are damped. Actually, 
there is a more dramatic effect. If all neurons have fired after a finite time t then all perturbed trajectories 
collapse onto the trajectory of V after t + 1 iterations. This loss of initial condition in a finite time is 
typical for IF models and is due to the reset of the membrane potential to a fixed value. For a discussion 
on IF model dynamics when this condition is relaxed see [113j . See also [79l HOlj . 

Initial conditions sensitivity. On the opposite, assume that there is a time, to, such that the image 
of the ball B(V, e) intersects S. By definition, this means that there exists a subset of neurons . . . , ik} 
and V £ B(V,e), such that Z(Vi(t )) ^ Z(V-(to)), i £ {ii, . . . , ik}. For example, some neuron does 
not fire when not perturbed but the application of an e-perturbation induces it to fire (possibly with 
a membrane potential strictly above the threshold). This requires obviously this neuron to be close 
enough to the threshold. Clearly, the evolution of the unperturbed and perturbed trajectory may then 
become drastically different. Indeed, even if only one neuron is lead to fire when perturbed, it may 
induce other neurons to fire at the next time step, etc . . . , inducing an avalanche phenomenon leading to 
unpredictability and initial condition sensitivity! 10 !. 

It is tempting to call this behaviour "chaos" , but there is an important difference with the usual notion 
of chaos in diffcrentiable systems. In the present case, due to the sharp condition defining the threshold, 
initial condition only occurs at sporadic instants, whenever some neuron is close enough to the threshold. 
Indeed, in certain periods of time the membrane potential typically is quite far below threshold, so that 
the neuron can fire only if it receives strong excitatory input over a short period of time. It shows then 
a behaviour that is robust against fluctuations. On the other hand, when membrane potential is close to 
the threshold a small perturbation may induce drastic change in the evolution. 

Stability with respect to small perturbations. Therefore, depending on parameters such as the 
synaptic efficacy, the external current, it may happen that, in the stationary regime, the typical tra- 
jectories stay away from the singularity set, say within a distance larger than e > 0. Thus, a small 
perturbation (smaller than e) does not produce any change in the evolution. At a computational level, 
this robustness leads to stable input-output transformations. 

On the other hand, if the distance between the set where the asymptotic dynamics liveiO and the 
singularity set is zero (or practically, very small) then the dynamics exhibit initial conditions sensitivity, 
and chaos. Typically a measure of this "distance" is given by |38j : 

9 Since time is discrete a neuron can fire and nevertheless satisfy this condition. 

10 This effect is well known in the context of synfire chains 2 3 4 91 or self-organized criticality |21| . 

11 Namely, the tj-limit set, f2, which is the set of accumulation points of F^A^), where P' (M.) is the mapping defining 
the dynamics (eq. 1141 1'). Since M. is closed and invariant, we have f2 = r)t=o F£(.M). In dissipative systems (i.e. a volume 
element in the phase space is dynamically contracted), the oj-limit set typically contains the attractors of the system. 
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d(fl,S)= inf inf min \VAt) - 61 (44) 
ven 4>o i=i...N 



where Q, is the w-limit set. 



Generic dynamics. Now, the following theorem holds [45] . 
Theorem 1. Ifd(Q,S) > then 

1. is composed of finitely many periodic orbits with a finite period, 

2. There is a one-to-one correspondence between a trajectory on fl and its raster plot, 

3. There is a finite Markov partition. 

Note however that d(fl, S) > is a sufficient but not a necessary condition to have a periodic dynamics. 
The main role of the condition d(fl,S) > is to avoid situations where the membrane potential of some 
neuron accumulates on 8 from below (ghost orbits). This corresponds to a situation where the membrane 
potential of some "vicious" neuron fluctuates below the threshold, and approaches it arbitrary close, with 
no possible anticipation of its first firing time. This leads to an effective unpredictability in the network 
evolution, since when this neuron eventually fire, it may drastically change the dynamics of the other 
neurons, and therefore the observation of the past evolution does not allow one to anticipate what will be 
the future. In some sense, the system is in sort of a metastable state but it is not in a stationary state. 

Now, assuming that conductances depend on past time only via a finite time horizon, one can show 
that, 

Theorem 2. Generically, in a probabilistic and topological sense, d(Q,S) > 0. 
(see [35] for the proof). 



Discussion Though the previous results suggests that dynamics is rather trivial since the first item 
tells us that dynamics is periodic, periods can however be quite long, depending on parameters. Indeed, 
following [38j an estimate for an upper bound on the orbits period is given by: 

T ~ 2 N io E «^» (45) 

where < 7 > denotes the value of 7 averaged over time and initial conditions. Though this is only 
an upper bound this suggests that periods diverge when d(fl,S) — > 0. This is consistent with the fact 
that when d(fi, S) is close to dynamics "looks chaotic" . Therefore, d(£l, S) is what a physicist could 
call an "order parameter", quantifying somehow the dynamics complexity. The distance d(Q,S) can be 
numerically estimated as done in [38j [45] . 



Let us give an example of application of this result. Consider model ()41|) where the synaptic weights 
are drawn at random with a Gaussian distribution A/"(^, ^-), in the same spirit as in section [2~Tl We have 
sketched the average value d(£l,S), averaged over the distribution of the Wij's, as a function of a, when 
J = and 7 is fixed. The curve of d(tt,S), as a function of a, delimits 3 regions. Region 1 corresponds 
to "neural death" (all neurons stop firing after a finite time); region 11 to a regime indistinguishable from 
chaos where the period of orbits are quite larger than what can be measured numerically; region 111 is 
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E[d(£2,S)] 
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Figure 4: Sketch of the dynamical regimes exhibited by model (1411 when synaptic weights are drawn at random with 

a Gaussian distribution J\f(0, (drawn from |38|). The expectation of d(Q,S) under the Wij's distribution, E[<i(fi,S)] 
, is drawn. It defines three regions. Region I corresponds to "neural death" (all neurons stop firing after a finite time); 
region II to a regime indistinguishable from chaos where the period of orbits are quite larger than what can be measured 
numerically; region III is a region where periodic orbits can be numerically detected. 



a region where periodic orbits can be numerically detected. This transition is reminiscent of the one 
exhibited in [110] for an isolated neuron submitted to a periodic excitation, but the present analysis hold 
at the network level. 

Let us now discuss the second item of theorem [TJ It expresses that the raster plot is a symbolic coding 
for the membrane potential trajectory. In other words there is no loss of information about the dynamics 
when switching from the membrane potential description to the raster plot description. This is not true 
anymore if d(fl, S) = 0. This issue, as well as the existence of a Markov partition, is used in section [31 

2.3 Conclusion 

In this section we have shown two examples of classical neural networks models, where the use of com- 
bined techniques from dynamical systems theory, statistical physics and probability theory allows the 
characterization of the dynamical regimes generically occurring. Moreover, considering random and in- 
dependent synaptic weights Wij 's we have been able to obtain a phenomenological "bifurcation diagram" 
where one replaces the overwhelming number of control parameters (TV 2 synaptic weights plus additional 
parameters defining the external current) by a small set of statistical parameters controlling the proba- 
bility distribution of the W^'s (mean and variance). This diagram characterizes the average behaviour 
of many different copies of the neural network when the W^ 's are drawn at random with a specific value 
of their mean and variance. It does not tell us what will be the typical behaviour of a given network 
(i.e. a given realization of the Wij's). Moreover, for the mean-field approach reported in section |2~T1 the 
bifurcation map corresponds to taking the limit N — *■ oo where, e.g. the transition to chaos is easy to 
represent since it is sharp. The situation is radically different for finite N where the "edge of chaos" 
associated with the transition by quasi-periodicity is rather complex and results from the overlapping of 
Arnold tongues [1211 [73] . For the gIF model, theorem [T] and [2] hold for generic values of the synaptic 
weights Wij's hence they apply to the huge space of parameters 7. Moreover, they characterize generic 
behaviours both in a topological and probabilistic sense. However, to figure out how d(tt, S) looks like 
we focused actually on the same situation as in section [2 . 1 1 where the Wij's are drawn at random, inde- 
pendently, where we study the effect of their variance of the average value of d(fl,S). It seems possible 
to have an analytic expression of G?(f2,5), but this requires to take the "thermodynamic limit" N — > 00 
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(Cessac and Touboul, in preparation). 

Thus, it appears clearly that these approaches are limited 

1. By the assumption of independence of the W^ 's. 

2. By the necessity of taking the limit N — > oo to obtain analytic expression. 
These limitations are further discussed in the conclusion section [5J 

3 Spikes trains statistics 

As we have seen in section [TJ neuronal activity is manifested by emission of spike trains having a wide 
variety of forms (isolated spikes, periodic spiking, bursting, tonic spiking, tonic bursting, etc) [551 [5(31 ITTTj . 
depending on physiological parameters, but also on excitation coming either from other neurons or from 
external sources. From these evidences, it seems natural to consider spikes as "information quanta" or 
"bits" and to seek the information exchanged by neurons in the structure of spike trains. Doing this, 
one switches from the description of neurons in terms of membrane potential dynamics, to a description 
in terms of spikes trains and raster plots. Though this change of description raises many questions it is 
commonly admitted in the computational neuroscience community that spike trains contain the "neural 
code". 

Admitting this raises however other questions. How is "information" encoded in a spike train: rate 
coding [5], temporal coding |167j . rank coding [1401 [62] . correlation coding [105] ? How to measure the 
information content of a spike train ? There is a wide literature dealing with these questions [1361 11041 
IT41 11351 [5J 11601 [711 1138] , which are inherently related to the notion of statistical characterisations of spike 
trains, see |145l [501 ITS] and references therein for a a review. As a matter of fact, a prior to handle 
"information" in a spike train is the definition of a suitable probability distribution that matches the 
empirical averages obtained from measures. 

3.1 Spike responses of neurons 

Neurons respond to excitations or stimuli by finite sequences of spikes. Thus, the dynamical response R of 
a neuronal network to a stimuli S (which can be applied to several neurons in the network), is a sequence 
uj(t) . . . uj(t + t) of spiking patterns. "Reading the neural code" means that one seeks a correspondence 
between responses and stimuli. However, the spike response does not only depend on the stimulus, but 
also on the network dynamics and therefore fluctuates randomly. Thus, the spike response is sought 
as a conditional probability P(R\S) |145j and "reading the code" consists of inferring P(S\R) e.g. via 
Bayesian approaches, providing a loose dictionary where the observation of a fixed spikes sequences R 
does not provide a unique possible stimulus, but a set of stimuli, with different probabilities. Having 
models for conditional probabilities P(R\S) is therefore of central importance. For this, one needs a good 
notion of statistics. 

These statistics can be obtained in two different ways. Either one repeats a large number of experi- 
ments, submitting the system to the same stimulus S, and performs a sample averaging. This approach 
relies on the assumption that the system has the same statistical properties during the whole set of 
experiments (i.e. the system has not evolved, adapted or undergone bifurcations meanwhile). Or, one 
performs a time average. For example, to compute P(R\S), one counts the number of times n(R,T,uj) 
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when the finite sequence of spiking patterns R, appears in a spike train uj of length T, when the network 
is submitted to a stimulus S. Then, the probability P(R\S) is estimated by: 

P(R\S)=lun ^>1. 

T— >oo 1 

This approach implicitly assumes that the system is in a stationary state. 

The empirical approach is often "in-between" . One fixes a time window of length T to compute the 
time average and then performs an average over a finite number AT of experiments corresponding to 
selecting different initial conditions. In any case the implicit assumptions are essentially impossible to 
control in real (biological) experiments, and difficult to prove in models. So, they are basically used as 
"working" assumptions. To summarise, one observes, from AT repetitions of the same experiment, Af 
raster plots u> m , m — 1 . . . Af on a finite time horizon of length T. From this, one computes experimental 
averages allowing to estimate P(R\S) or, more generally, to estimate the average value, (0), of some 
prescribed observable 4>{lo). These averages are estimated by : 

^ T) = iEEfr)- ( 46 ) 

m=l t=l 

Typical examples of such observables are (f>(u>) = Wi(0) in which case (<fi) is the firing rate of neuron i; 
4>(uj) = uJi(0)u>j(Q) then ((f)) measures the probability of spike coincidence for neuron j and i; <f)(o->) — 
u)i(T)u)j(0) then ((/)) measures the probability of the event "neuron j fires and neuron i fires r time step 
later" (or sooner according to the sign of t). In the same way P(R\S) is the average of the indicatrix 
function X-r( w ) = 1 if w G i? and otherwise, the statistics being performed when the neuronal network 
is submitted to S. Note that in (|46|) we have used the shift <r* for the time evolution of the raster plot. 
This notation is more compact and more adapted to the next developments than the classical formula, 
reading, e.g., for firing rates £m=i E*Li 0( w m(*))- 

This estimation depends on T and Af. However, one expects that, as Af, T — > oo, the empirical 
average (j)^' 7 "* converges to the theoretical average (0), as stated e.g. from the law of large numbers. 
Unfortunately, one usually does not have access to these limits, and one is lead to extrapolate theoretical 
averages from empirical estimations. The main difficulty is that these observed raster plots are produced 
by an underlying dynamics which is usually not explicitly known (as it is the case in experiments) or 
impossible to fully characterise (as it is the case in most large dimensional neural networks models). 
Thus, one is constrained to propose ad hoc statistical models. As a matter of fact, the choice of a statis- 
tical model always relies on assumptions. Here we make an attempt to formulate these assumptions in a 
compact way with the widest range of application. These assumptions are compatible with the statistical 
models commonly used in the literature like Poisson models or Ising like models a la Schneidman and 
collaborators [155] . but lead also us to propose more general forms of statistics. Moreover, our approach 
incorporates additional elements such as the consideration of neurons dynamics, and the fact that this 
dynamics severely constrain the set of admissible raster plots, £ 7 . This last issue is, according to us, 
fundamental, and, to the best of our knowledge, has never been considered before in this field. 

On this basis we propose the following definition. Fix a set (j>i, I — 1...K, of observables, i.e. 
functions £ 7 — > IR which associate real numbers to sequences of spiking patterns. Assume that the 

empirical average ([46]) of these functions has been computed, for a finite T and Af, and that <f>i = C). 
A statistical model is a probability distribution v on the set of raster plots such that: 
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1. — 1, i.e. the set of non admissible raster plots has a zero ^-probability. 

2. v is ergodic for the left-shift <t 7 . 

3. For all I = 1 . . . K, v(4>i) — Ci, i.e., v is compatible with the empirical averages. 

Note that item 2 amounts to assuming that statistics are invariant under time translation. On 
practical grounds, this hypothesis can be relaxed using sliding time windows. This issue is discussed in 
more details in [40 . Note also that v depends on the parameters 7. Assuming that v is ergodic has the 
advantage that one does not have to average both over experiments and time. It is sufficient to focus on 
time average for a single raster plot, via the time-empirical average: 



3.2 Raster plots statistics. 

A canonical way to construct statistical models comes from statistical physics |103j . This approach has 
been introduced for spike train analysis by [155] and generalised in [40j . According to item (l)-(3) we are 
seeking a probability distribution v which matches the constraints v{4>i) = Cj, I = 1 . . . K, where v(<j>i) 
is the average of </>; under v. We want to stick on these constraints, imposed by experimental results, 
without adding any other hypothesis. In the realm of statistical physics this amounts to maximising 
the statistical entropy under the constraints v{4>i) = Cj 5 I = 1 ■ ■ . K. In the context of the so-called 
thermodynamic formalism of ergodic theory, which is a quite powerful tool to handle such statistical 
problems, this amounts to solving the following variational principle: 

P[ip}= sup (h [v] + v [V]), (48) 

where the set of invariant (stationary) measures for the dynamics and h is the entropy rate. We 

have introduced a "potential" , 

K 

v> = 5>^, ( 49 ) 

1=1 

where the Az's are adjustable Lagrange multipliers. A probability measure which realises the supre- 
mum, i.e. 

P[ip] = h [vj,] + , 

is called an "equilibrium state" . The function P [if>] is called the "topological pressure" in the realm of 
ergodic theory, and "thermodynamic potential" (free energy, free enthalpy, pressure) in statistical physics. 
Note that ergodic theory imposes less constraints on dynamics than statistical physics (the microscopic 
dynamics does not need to be Hamiltonian) . From the topological pressure one computes the moments 
of the distribution v^. In particulaJ^l. 

^ - M<h). (50) 



2 This relations assumes that P [tp] is differentiable, i.e. that the system is away from a phase transition. 
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This relation fixes the value of the Lagrange multipliers A/ in order to have (<fri ) = C; . 

Moreover, in "good cases" (e.g. uniformly hyperbolic dynamical systems), equilibrium states are also 
Gibbs states [231 [2^1 [TT^l 11391 H7] , A Gibbs state, or Gibbs measure, is a probability measure such that, 
one can find some constants ci, ci with < c\ < 1 < ci such that for all n > 1 and for all to: 



ci < 



M0 (uj G H „_ 1 ) 

exp(-nP [</>] +S'(")^(w)) 



<C 2 , 



(51) 



where S^ip(uj) = ^2™ = q ^{u^oj) and where we denote by [uj] n _ 1 a cylinder set of length n, namely the 
set of raster plots io' such that w'(i) = u>(t), t = . . . n — 1. Basically, this means that the probability 

that a raster plot starts with the bloc [w] 0n _j behaves like cxp ^ s ^ ) ^^"- ) - ) . One recognises the classical 
Gibbs form where space translation in lattice system is replaced by time translation (shift <j 1 ^) and where 
the normalisation factor Z n is the partition function. Note that P = limsupj^^ ilogZ„, so that 
P [i/>] is indeed the formal analog of a thermodynamic potential (like free energy) . 

In this context, the probability of a spiking pattern block R — [u;] n l of length n corresponding to 
the response R to a stimuli S "behaves like" (in the sense of eq. ([ST]) ): 



P[R\S] =v[ujE R\S] 



Z n [Ai(S),...,A,(5)] 



exp 



K 



5>(S)5>(^- 



.1=1 



t=0 



(52) 



where the Az's depend on the stimulus S. Obviously, for two different stimuli the probability P(R\S) may 
drastically change. 



3.3 Examples. 



(T) 

Firing rates. If 4>i{uS) — u>i(0), then 7r^ {4>i) = r\ is the average firing rate of neuron I within the time 
period T. Then, the corresponding statistical model is a Bernoulli distribution where neuron I has a 
probability ri to fire at a given time. The probability that neuron I fires k times within a time delay n 
is a binomial distribution and the inter-spike interval is Poisson distributed [77] . 

Spikes coincidence. If cf)i(u>) = <j>(i j)(w) = Wj(0) Wj(0) where, here, the index I is an enumeration for 
all (non-ordered) pairs then the corresponding statistical models has the form of an Ising model, 

as discussed by Schneidman and collaborators in [155, 170 . As shown by these authors in experiments 
on the salamander retina, the probability of spike blocs estimated from the "Ising" statistical model fits 
quite better to empirical date than the classical Poisson model. 



Enlarged spikes coincidence. As a generalisation one may consider the probability of co-occurrence 
of spikes from neuron i and j within some time interval r. The corresponding functions are 4>i(u>) = 
Wi(0)cjj-(r) and the probability of a spike bloc R reads: 



P[R\S] 



1 



Zn [Ai,i(5), . 

An example is provided in section f5. 31 



*N,N 



(S)] 



exp 



X>«(S)E^(*W* + t) 

i<j t=0 
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Further generalisations can be considered as well. 



Generalised Integrate and Fire models Due to their particular structure and especially the fact that 
generically a Markov partition exists, gIF models of type (|4"Tj) are explicit examples where this theory 
gives striking results (see [JU] for details and section [5] for an application to the effect of synaptic plasticity 
to spike trains statistics.) 



3.4 Validating a statistical model 

There are currently huge debates on the way how brain encodes information. Are frequency rates sufficient 
to characterise the neural code [175] ? Are pair correlations significant ? Do higher order statistics matter 
? Actually it might be that the answer depend on the brain process under consideration and some peoples 
actually believe that "brain speaks several languages and speak all of them at the same time" (Franck 
Grammont, private communication. For a nice illustration of this see |80j ) . These questions are inherently 
linked to the notion of (i) finding statistical models; (ii) discriminate several statistical models and select 
the "best one". 

Let us consider an illustrative example, i.e. the question: are correlations significant ? Answering this 
question is a crucial issue for biologists/experimentalists [156} 11411 H42j . Note that it has absolutely no 
meaning to try and answer this question from empirical data when considering "the brain" as a whole. 
But, as emphasised by [148J, there is maybe some hope to make one step forward when considering small 
neural assemblies (e.g. small pieces of retina). 

Moreover this question has no "absolute" answer but a relative answer in the following sense. Let us 
consider the 1st order potential: 

N 



thus only taking firing-rates into account, "against" the 2nd order potential: 

N N T B 

■02 M = Y Xi LJj (°) + Y Y Ai J' T U> i(°) UJ j( T )> 
i=l i,j=l t=-T s 

where T s is a characteristic time scale. This potential form takes both firing-rate and correlations into 
account. 

The realm of thermodynamic formalism offers a numerically tractable way to compare the statistical 
models related to these two potentials. The relative entropy or Kullback-Leibler divergenc^H between a 
Gibbs measure and a stationary measure /x is given by [1121 l4"r?I |4"T] : 



h [n\v^) = P [-0] - j i/)d(i - h(jj). 



13 Let be two invariant measures both denned on the same set of admissible raster plot S T . The relative entropy (or 
Kullack-Leibler divergence) between fi and v is: 



h{p\v) = limaup- Yl t* (Mo,b-i) 1o S 
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In the present case, we are given an empirical measure, ttIj , (see eq. ([47]) ). obtained from experiments. 
To discriminate between the two potentials xp 1 ,if> 2 a possible criterion consists of choosing the potential 
which minimises the relative entropy of the corresponding Gibbs measure with respect to the empirical 
measure. Namely, if there is To > such that, for all T > To : 

M4 T) I^)<M^ T) I^ 2 ) (54) 

then ip 1 is considered as a better statistical model than if> 2 . The nice thing is that, using the thermody- 
namic formalism, one can develop algorithms allowing such a comparison (Cessac, Vasquez, Vieville, in 
preparation) . 

3.5 Conclusion 

When analysing spike train statistics, one is lead to propose several statistical models corresponding to 
distinct hypothesis. For example, characterising inter-spike interval distribution by a homogeneous Pois- 
son process ultimately corresponds to assuming that correlations between neurons and time correlations 
are irrelevant and that frequency rates are sufficient to characterise statistics. In our presentation this 
amounts to considering a Gibbs potential of form J2iLi AjWi(O). More general forms can be proposed as 
well. But this leads to two fundamental questions: 

1. How to discriminate statistical models from empirical data ? This is a crucial issue, whose tractabil- 
ity was deeply raised by and Roudy and his collaborators in a very recent paper |148j . Many criteria 
used in the literature rely on the computation of the Kullback-Leibler divergence. We have shown 
how the use of the thermodynamic formalism, relying on a safe recipe from statistical mechanics, 
could allow to compute this quantity from data. 

2. Instead of defining the statistical model from ad hoc observables, is it possible to propose a canonical 
form relying on some generic principle ? This issue is addressed in section 15.31 where we consider 
the effect of synaptic plasticity on spike trains statistics. 

4 Interplay between synaptic graph structure and neurons dy- 
namics. 

4.1 Causal actions 

Since synapses are used to transmit neural fluxes (spikes) from a neuron to another one, the existence 
of synapses between a neuron (A) and another one (B) is implicitly attached to a notion of "influence" 
or causal and directed action. However, as we saw, a neural network is a highly dynamical object and 
its behavior is the result of a complex interplay between the neurons dynamics and the synaptic network 
structure. Moreover, the neuron B receives usually synapses from many other neurons, each them being 
"influenced" by many other neurons, possibly acting on A, etc... Thus the actual "influence" or action of 
A on B has to be considered dynamically and in a global sense, by considering A and B not as isolated 
objects, but, instead, as entities embedded in a system with a complex interwoven dynamical evolution. 
In this context it is easy to imagine examples where there is a synapse from A to B but no clear cut 
influence, or, in the opposite, no synapse and nevertheless an effective action. 
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Thus, one has to consider topological aspects (such as the feedback circuits) and dynamical aspects. 
One way of doing this is to compute correlations between neurons (cross-correllogramms). However, 
correlations functions do not provide causal information. A strong correlation between A and B at time 
t does not tell us if A acts on B or if B acts on A (note in particular that Cab(£) = Cba(— *))■ 

Another way consists of exciting neuron A, say with a weak signal, and observe the effects on £>, 
e.g. by comparing its evolution with and without the signal applied on A. A natural choice for an 
excitatory signal is a periodic signal, with a tunable frequency. Thus, the response function, drawn 
versus frequency, provides similar information as the complex susceptibility in Physics. In particular, 
peaks in the susceptibility corresponds to resonances, that is, a response of maximal amplitude. These 
resonances can be used to provide an effective, frequency dependent notion of network structure, as we 
now show. 

4.2 A simple but non trivial example 

Consider the model (|32|) where neurons are represented by frequency rates. As we saw in section [2.1.3l this 
model exhibit, in finite dimension N a generic transition to chaos by quasi-periodicity, when increasing 
the non-linearity of the sigmoidal transfer function S. 

Signal propagation and effects of non-linearity. Assume that this system is in the chaotic regime. 
Note that the corresponding Fourier spectrum is not flat but contains peaks (resonances) reminiscent of 
the transition by quasi-periodicity [13] . Assume now that we superimpose upon the membrane potential 
Vj(t) of the neuron j a small external signal £j(t). Does this signal have an effect which propagates inside 
the network ? and how ? Because of the sigmoidal shape of the transfer functions the answer depends 
crucially, not only on the connectivity of the network, but also on the value of the T4's. Assume, for the 
moment and for simplicity, that the time-dependent signal £j (t) has variations substantially faster than 
the variations of Vj. Consider then the cases depicted in Fig. \5[ In the first case (a) the signal £j(t) is 
amplified by S, without distortion if is weak enough. In the second case (Fig. [SJd), it is damped 
and distorted by the saturation of the sigmoid. More generally, when considering the propagation of 
this signal from the node j to some node i one has to take into account the level of saturation of the 
nodes encountered in the path, but the analysis is complicated by the fact that the nodes have their own 
dynamical evolution (Fig. Hb). 

Tangent space splitting. In this context we would like to measure the average "influence" of neuron 
A on neuron B (namely how a weak signal applied on A perturbs on average the state of B) 1 including the 
effects of the nonlinear dynamics. Typically, in dissipative systems, such as neural networks models, where 
volume in the phase space is dynamically contracted, dynamics asymptotically settle "onto" attractorJ^l. 
Examples of attractors are stable fixed points, or stable periodic orbits. Chaotic attractors have moreover 
the following property. While dynamics transverse to the attractor is contracting (corresponding precisely 
to the attractivity property), dynamics "parallel" to the attractor is expanding, corresponding to initial 
condition sensitivity. In other words, the tangent space of attractor points can be split into a contracting 
and an expanding part (see fig. |BJ). 

14 Let JVl be the (compact) phase space of the dynamical system. A set A £ M is called an attractor if it is invariant 
(F t (A) = A) and if if there exists an open set U G -M such that A = n t >oF*(W). This definition affords several non 
equivalent extensions and variants 11801 IT291 [65] [54] [55] [1091 
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Figure 5: Nonlinear effects induced by a transfer function with a sigmoidal shape on signal transmission. Fig. [5^ (left). 
Amplification. Fig. [5p (middle). Saturation. Fig. 0; (right). The propagation of a signal along a path in the network 
depends not only on the weights of the links but also on the level of saturation of the nodes that the signal meets. The 
level of saturation depends on the current state of the node (schematically represented as a red point in the figure). This 
state evolves with time. 



Structural properties of chaotic attractors are usually characterized by statistical quantities such as 
Lyapunov exponents. There is indeed a natural notion of average in chaotic systems related to the Sinai- 
Ruelle-Bowen measure p (SRB) |159| l25l 11501 which is obtained as the (weak) limit of the Lebesgue 
measure p, under the dynamical evolution 1 ^: 

p = lim F"/i. 

n — >+oc ' 

where F™ is the n-th iterate of F 7 . Following an orbit upon the attractor it is possible to characterize 
the average expansion and contraction rates for this orbit via Lyapunov exponents. A positive Lyapunov 
exponent indicates local expansion while a negative one indicates contraction. 




Figure 6: Sketch of an attractor and of the decomposition of the tangent space into contracting and expanding directions. 
A perturbation in the contracting space leads to a trajectory which converges back to the attractor. A perturbation in 
the transverse direction leads to a trajectory staying onto the attractor but which separates from the mother trajectory. 
This representation is made for a discrete time system. Note that the regular spacing between points is an artefact of the 
representation. 



In the following we will assume that all Lyapunov exponents are bounded away from zeroF^I. Then for 
each V € supp p, where supp p is the support of p, there exists a splitting © such that E^ , the 

15 A crucial property is that a SRB measure has a density along the unstable manifolds, but it is singular in the directions 
transverse to the attractor. 

16 This formalism requires, on rigorous grounds, that the system is uniformly hyperbolic 11511 . and examples of diverging 



35 



unstable space, is locally tangent to the attractor (the local unstable manifold) and E-^ , the stable space, 
is transverse to the attractor (locally tangent to the local stable manifold). Let us emphasize that the 
stable and unstable spaces depend on V (while the Lyapunov exponents are p almost surely constant). 

Let us consider a point V upon the attractor and make a small perturbation 8-y. This perturbation 
can be decomposed as #v = + wnere ^ -^V an< ^ ^v"* e -^V ■ * s l° ca Hy amplified with an 
exponential rate (given by the largest positive Lyapunov exponent). On the other hand 6-^ is damped 
with an exponential speed (given by the smallest negative Lyapunov exponent) 

Linear response. Assume now that we superimpose a signal of weak amplitude, considered as an 
external current, upon some neuron (fc) in such a way that the dynamics is still chaotic (with only a tiny 
variation of the Lyapunov exponents). For simplicity, we suppose that the signal does not depend on 
the state of the system, but we can consider this generalization without difficulty (linear response still 
applies in this case, but equations (|55|) . ([56]) do not hold anymore). Denote by £ the vectoi0 
The new dynamical system is described by the equation: 



The weak signal £(t) may be viewed as a small perturbation of the trajectories of the unperturbed system 
([32]) . At each time this perturbation has a decomposition £(t) = £^(f) + ^ u \i) on the local stable and 
unstable spaces. The stable component is exponentially damped. The unstable one ^ u '(t) is 

amplified by the dynamics and then scrambled by the nonlinear terms. Consequently, it is impossible to 
predict the long term effect of signal £(t) on the global dynamics. 

This is true for individual trajectories. However, the situation is substantially different if one considers 
the average effect of the signal, the average being performed with respect to the SRB measure p of the 
unperturbed system. Indeed, as an application of the general theory [151] . it has been established in 
[42], l43l [44] that the average variation 8y i it) of the membrane potential Vi under the influence of the 
signal is given, to the linear order, by: 



We used the shortened notation < > for the average with respect to p. In this expression Xij (c) ar e the 
matrix elements of : 



Thus x( tJ ) is a matrix representing the average value of the iterate a of the Jacobian matrix. Let us 
note that the fact that xi a ) remains bounded for a — ► oo is not a trivial result because DF^ diverges 
exponentially with a. The convergence of xi a ) nas been rigorously shown by Rucllc under the hypothesis 
of uniform hyperbolicity. It results from the exponential correlation decay (mixing) in the unstable 
directions and on the exponential contraction. 

susceptibility can be exhibited for the logistic map |152| or the Henon map [37| . Also, it does not hold at a bifurcation 
point where susceptibility can diverge. On practical grounds we require that all Lyapunov exponents are bounded away 
from zero. 

17 Though this term acts in equation H32H as an external current, we use the notation £ throughout the paper, to 
distinguish between an arbitrary external current (l( e:Et )) and some specific stimulus intended to carry some "information" 
in the network. 



y(* + i) = F 7 v(t) +€(t). 



oo 




(55) 
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This means that, provided that £(t) is sufficiently small, and for any smooth observable A, the 
variation < A > t — < A > is proportional to up to small nonlinear corrections. In other words, pt 
is differentiable with respect to the perturbation. The derivative is called the linear response. 

Causal circuits. In the case of the dynamical system (plilj) one can decompose Xij( T ) as : 



The sum holds on each possible paths Jij(r), of length t, connecting the neuron fco = j to the neuron 
k r = i, in t steps. One remarks that each path is weighted by the product of a topological contribution 
depending only on the weight W%j and a dynamical contribution. Since, in the kind of systems we 
consider, functions S are sigmoid, the weight of a path 7^ (r) depends crucially on the state of saturation 
of the neurons k , . . . , k T -i at times 0, . . . , r — 1. Especially, if S'(Vk l _ 1 (I - 1)) > 1 a signal is amplified 
while it is damped if S'(Vk l _ 1 {I — 1)) < 1. Thus, though a signal has many possibilities for going from j 
to i in t time steps, some paths may be "better" than some others, in the sense that their contribution to 
Xij( T ) is higher. Therefore eq. (|F7)) , which quantifies the intuition raised in fig. [5J underlines a key point. 
The analysis of signal transmission in a coupled network of dynamical neurons with nonlinear transfer 
functions requires to consider both the topology of the interaction graph and the nonlinear dynamical 
regime of the system. 

As a remark note that since the derivatives S' in (|57[) are bounded by some constant (proportional to 
g), one can bound the Jacobian matrix component |-D-Fy| by some Xij. This provides a bound on (|57[) 
which resembles very much to an expression obtained by Afraimovich and Bunimovich in [6j (lemma 3) 
from which they derive, using the thermodynamic formalism, a topological pressure characterizing the 
stability of the dynamical system. Actually, their analysis fully applies here when the attractor is a fixed 
point and their theorem 1 typically provides a parametric condition for the stability of the fixed point. 
Note that then eq. (|57| reduces to Xij{ T ) = S 7i -( T ) Y\i=i Wfcjfcj-i an d expresses the Jacobian matrix DF T 
at the fixed point, in terms of graph loops. This can be related, to the so-called cyclic expansions used in 
dynamical system and ergodic theory (see http://chaosbook.org/ for a very nice presentation). Actually 
we believe that Afraimovich and Bunimovich approach can be extended to our case also in the case of 
chaotic dynamics. But one has to use a double cyclic expansion: on the loops of the graph, and on the 
unstable periodic orbits which can be used to approximate the SRB measure (Cessac, in preparation). 

Complex susceptibility. The existence of this linear response theory opens up the way to applications 
involving chaotic neural networks used as a linear filter. Indeed eq. (|55[) describes a linear system which 

transforms an input signal £(t) of small amplitude into an output signal (Vi{t) — Vi(t) \ according to a 

standard convolution product. In particular, if the external signal is chosen as: 



(where e^ is the unit vector in direction j), then the response of the system is also harmonic with : 




(57) 



ce 
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Vi{t) - Vi{t) 
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where the frequency-dependent amplitude: 

oo 

*«M = 5>ii(*)e fa " r (58) 

(7=0 

is called the complex susceptibility. In ref. (42] a method have been conceived and implemented allowing 
to compute \ij ( w ) numerically. The knowledge of the susceptibility matrix is very useful as it enables 
one to detect resonances, i.e. frequencies for which the amplitude response of the system to a periodic 
input signal is maximum. In fact the existence of a linear response implies that Xij(w) is bounded for 
all to £ [0,27r]. Moreover, in view of eq. (|58p . it is analytic in the complex upper plane. On the other 
hand, Xjj(w) can have poles within a strip in the lower half plane, e.g. in loq — iA, A > 0. In this case, 
and if A is small, the amplitude 1x^(^)1 exhibits a peak of width A and height IXijC^o)! which can be 
interpreted in the present context as follows: when unit j (whose state varies chaotically due to the global 
dynamics) is subjected to a small periodic excitation at frequency ujq and amplitude e then the average 
response of unit i behaves periodically with same frequency and amplitude e|xij(wo)| which is maximal 
in a frequency interval centered about too. 

An example of resonances The following case has been analysed in [33] for details. This is a sparse 
network where each unit receives connection from exactly K = 4 other units. The corresponding network 
is drawn in Fig. [JJi. Blue stars correspond to inhibitory links and red crosses to excitatory links. In this 
example the unit 7 is a "hub" in the sense that it sends links to most units, while 0, 2, 3 or 5 send at 
most two links. 

In figure Et), we have represented the modulus of the susceptibilities for all pairs and different 
frequencies u). This provides a notion of "causal connectivity", related to linear response, which departs 
strongly from the connectivity provided by weights matrix. 

Computing the susceptibility one obtains the curves shown in FigUJ;. Some resonance peaks are rather 
high (~ 20) corresponding to an efficient amplification of a signal with suitable frequency. It is also clear 
that the intensity of the resonance has no direct connection with the intensity or the sign of the coupling 
and is mainly due to nonlinear effects. For example, there is no direct connection from to 3 or 5 but 
nevertheless these units react strongly to a suitable signal injected at unit 0. Let us now compare the 
Fourier transform of the correlations function Cij(t) for the same pairs (Fig. [7Jd). One remarks that 
these functions exhibit less resonance peaks. This is explained in the context of Ruelle's linear response 
theory and is related to the decomposition of the linear response into stable and unstable contributions, 
related to the local splitting of the tangent space (see (42] [43] for more details). 

4.3 Conclusion. 

The previous analysis leads then us to propose a notion of "effective", frequency dependent, connectivity 
based on susceptibility curves. For a given frequency to, we plot the modulus of the susceptibility IXij^)! 
with a representation assigning to each pair i,j a circle whose size is proportional to the modulus. We 
clearly see in figure [7] that changing the frequency changes the effective network. 

All these effects are due to a combination of topology and dynamics and they cannot be read in 
the connectivity matrix W. Therefore, the example of neural networks treated in this section shows 
convincingly that the analysis of neural circuits requires a careful investigation of the combined effects 
induced by non-linear dynamics and topology of the synaptic graph. It also shows that the analysis of 
correlations provides less information than a linear response analysis. This is particularly clear when 
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looking at the resonances curves displayed by linear response and correlations functions. As we saw, this 
difference is well understood on theoretical grounds and has deep relations with salient characteristics 
of the nonlinear dynamics (saturation in the transfer function closely related to the refractory period) . 
Using linear response in neural networks is not new (see for example [145] and references therein), but 
the point of view adopted in the present section, is, we believe, less known and raises new interesting 
questions. 

For example, one may wander what would bring this approach in spiking neural networks, where the 
causal action from a neuron to another can be somewhat "directly" read in the timing of pre- and post- 
synaptic neurons spikes. Another remaining question is what would the use of linear response analysis 
tell us in neural networks having synaptic plasticity. This issue is briefly addressed in the next section. 

5 Dynamical effects of synaptic plasticity 
5.1 General context 

The notions of neural code and information cannot be separated from the capacity that neuronal net- 
works have to evolve and adapt by plasticity mechanisms, and especially synaptic plasticity. Therefore, 
understanding the effects of synaptic plasticity on neurons dynamics is a crucial challenge. Especially, 
addressing the effect of synaptic plasticity in neural networks where dynamics is emerging from collective 
effects and where spikes statistics are constrained by this dynamics seems to be of central importance. 
This issue is subject to two main difficulties. On one hand, one must identity the generic dynamical 
regimes displayed by a neural network model for different choices of parameters (including synaptic 
weights) . Some examples of such analysis have been given in section [2] On the other hand, one must 
analyse the effects of varying synaptic weights when applying plasticity rules. This requires to handle 
a complex interwoven evolution where neurons dynamics depends on synapses and synapses evolution 
depends on neuron dynamics. 

Effects of synaptic adaptation Three main classes of effects can be anticipated [59 j 11611 1162] . 

1. Structural effects. There is a first, evident, effect of synaptic plasticity: a rewiring of the neural 
network. However, this rewiring is not some random process where edges would be selected or 
removed independently of the history. Instead, it results from a complex process where edges are 
potentiated or depressed according to the neuron dynamics, which is itself depending on synaptic 
weights and external stimuli. The question is therefore whether one can nevertheless extract some 
general characteristics of the network structure evolution and what is the impact of this structure 
evolution on neural network behaviour. 

2. Dynamical effects. Changing the synaptic weights, which are parameters of the dynamical sys- 
tems (|13j) and (|14| . will obviously have an incidence on dynamics. These effects can be smooth or 
sharp (bifurcations). More generally, one expects period of smooth changes interrupted by sharp 
transitions (see fig. [5] for an example of this). Thus, adaptation drives the dynamical system along 
a definite path in the space parameters, which integrates the whole past, via synaptic changes. In 
this respect, we address a very untypical and complex type of dynamical system. This induces rich 
properties such as a wide synaptic adaptation-induced variability in the network response to a given 
stimulus, with the same set of initial synaptic weights, simply by changing the initial conditions. 
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3. Functional effects. This evolution typically arises when the system is submitted to inputs or 
stimuli which constrain neuron dynamics and thus synaptic evolution. For simplicity, let us think 
of synaptic plasticity in the restricted context of "learning" some input. Learning should result in 
the acquisition of a new ability. The network after learning should be able to "recognise" a learnt 
input, while this was not necessarily the case before learning. In this context, recognition can be 
manifested by a drastic change in the dynamics whenever a learnt input is presented. Moreover, 
this effect must be robust and selective. Examples of this are presented now. 

Coupled dynamics As an illustration we consider now the following coupled dynamics. Neurons 
are evolving according to (|14[) (we focus here on discrete time dynamics). We consider slow synapses 
dynamics. Namely, synaptic weights are constant for T consecutive dynamics steps, where T is large. 
This defines an "adaptation epoch" . At the end of the adaptation epoch, synaptic weights are updated 
according to @. This has the consequence of modifying neurons dynamics and possibly spike trains. 
The weights are then updated and a new adaptation epoch begins. We denote by t the update index 
of neuron states (neuron dynamics) inside an adaptation epoch, while r indicates the update index of 
synaptic weights (synaptic plasticity). Call X( r ) (t) the state of the neurons at time t within the adaptation 
epoch t (we use here the notation X instead of V since eq. (I5"9l holds for model having possibly more 
variables than the membrane potential for the definition of the neuron state). Let be the synaptic 

weights from neuron j to neuron at i in the r-th adaptation epoch. At the end of each adaptation epoch, 
the neuron dynamics indexes are reset, and X^ T+1 \o) — x[ T \T),i = 1...N. The coupled dynamics 



Recall that 7 = {W,I (ext) ) (see section OH) and 7^ T - ) is the set of parameters at adaptation epoch 
r. In the present setting the external current is used as a time constant external stimulus. We write it £ 
(see footnote [17]) . 

Let us discuss a few examples. 

5.2 Hebbian learning 

This example has been considered in [1611 1162] . The goal is to study the role of synaptic plasticity 
mechanisms inspired by Hebb's work 9Q\ and its generalisation, in a situation where the neural network 
is submitted to some specific stimulus over a long time. The main is to study the conjugated effects 
of stimulus action and synaptic plasticity mechanisms on the neural network dynamics. Specifically, we 
want to investigate whether these conjugated effects can lead the system to a state where it acquires 
some ability to "recognise" this stimulus. In the example developed below, this corresponds to drive the 
system, via an Hebb's- inspired modification of the synaptic weights, in a region in the parameters space 
where presenting the stimulus induces a bifurcation in the dynamics whereas this bifurcation didn't occur 
before synaptic adaptation. 

Let us before briefly explain what we mean by "Hebbian learning". D. Hebb has proposed in [90] 
a theory of behaviour based on the physiology of the nervous system. The most important concept to 
emerge from Hebb's work was his formal statement (known as Hebb's rule) of how learning could occur. 



writes: 




(59) 
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When an axon of cell A is near enough to excite a cell B and repeatedly or persistently takes part in 
firing it, some growth process or metabolic change takes place in one or both cells such that A 's efficiency, 
as one of the cells firing B, is increased. 

Many "learning rules" in neural networks are based on Hebb's observations plus a few well established 
facts. They rely upon a few recipes that can summarised as [57] : 

• Learning results from modifying synaptic connections between neurons. 

• Learning is local i.e. the synaptic modification depends only upon the pre- and post- synaptic 
neurons activity and does not depend upon the activity of the other neurons. 

• The modification of synapses is slow compared with characteristic times of neuron dynamics. 

• If either pre- or post- synaptic neurons or both are silent then no synaptic change takes place except 
for (exponential) decay which corresponds to forgetting. 

The first item implies that learning results in a modification of the Wij's. The second one basically 
says that the synaptic modification of Wij writes W[j — eh(W? ,mj,mi) where WL is the value of the 
synapses j — ► i after the learning rule has been applied. The numbers TOj (to.,) denotes the "state" 
or "activity" of the neuron i (j). How this "state" is defined vary according to the model. The third 
item implies then that e is a small parameter, whose inverse corresponds to the characteristic time for a 
significant change of Wij . If one assumes that h is a smooth function then one may simply consider a 
Taylor expansion of a generic regular function h. This gives, up to the second order in mi,mj. 

W-j = e (a oo + aiooWy + aoiom? + aooimj + aoiimjTOj + h.o.t.) 

where h.o.t. means "higher order terms" such as WijiniWij, etc.... Then, the fourth item implies aooo = 
and leads to introduce a parameter A = eaioo G [0, 1] that models passive "forgetting": if a synapse is 
not solicited its intensity decreases with a decay rate j. 



5.2.1 Coupled dynamics. 

On these bases, consider the model (|52")l where synaptic weights evolve according tj^l: 

N 

(60) 

where a is a small number, controlling the rate of synaptic plasticity. Here, one associates to the history 

(t) 

of neuron i s rate an activity index : 

™| T) = ^X>| T) M-^) 

where (t) is the firing rate of neuron i at time t in the adaptation epoch r, di € [0, 1] is a threshold. 

The neuron is considered active during synaptic adaptation epoch r whenever m\ T ^ > 0, and silent 
otherwise. Finally, H(x) is the Heaviside function. 

18 For another implementation of Hebbian rule see eq. (|10fl . 



42 



Let us now interpret this equation. The first term is conform to the recipes introduced in the previous 
section. The second term is positive whenever neuron j and i are active, increasing the synapse efficacy 
Wij (i.e. render it more positive if is excitatory and less negative if it is inhibitory). If j is active and 
i inactive the synapse efficacy decreases. Finally, if j is inactive the second term is zero. We emphasise 
that this one possible implementation among many others. 

Equations ([32]) & (|60|) define a dynamical system where two distinct processes (neuron dynamics 
and synaptic network evolution) interact with distinct time scales. This results in a complex interwoven 
evolution where neuronal dynamics depends on the synaptic structure and synapses evolve according to 
neuron activity. On general grounds, this process has a memory that is a priori infinite and the state of 
the neural network depends on the past history. 

5.2.2 Observed effects of Hebbian synaptic plasticity 

The effect of this Hebbian synaptic adaptation has been explored numerically in [59] and mathematically 
in [162j . Assume that the initial synaptic weights are chosen independently, at random, with a Gaussian 
distribution of mean ^ and variance jt- . This choice mimics a situation where no structure is imposed a 
priori in the correlations between synaptic weights. The idea is to see how synaptic adaptation changes 
this situation. Then, according to section (|2.1.3|) dynamics is chaotic provided the gain g of the sigmoid 
S is large enough. Starting from such a chaotic dynamics, the following effects have been observed, in 
correspondence with the three effects anticipated in the beginning of this section. 

Structural effects. The rewiring of the network by the synaptic adaptation rule (pO)) reveals a variation 
in the synaptic weights distribution, and an increase in weights correlations. This effect is evident but 
does not give significant hints to interpret the dynamical and functional effects described below )162j . 
Also, a computation of standard indicators in complex graphs analysis, when applied to the synaptic 
weights matrix, does not show any salient effect. The only important hint provided b y s ynaptic weights 
matrix analysis is an increase in the number and weights of positive feedback loops 19 !, which renders 
the system more cooperative, with a strong impact on dynamical [93] . As suggested in the section [4] 
the analysis of synaptic weights matrix is indeed not expected to provide deep insights on dynamical 
effects. On the opposite the analysis of Jacobian matrices reveals important properties. This is not 
surprising. A standard procedure for the analysis of nonlinear dynamical systems starts with a linear 
analysis. This holds e.g. for stability and bifurcation analysis but also for the computation of indicators 
such as Lyapunov exponents. The key object for this analysis are Jacobian matrices. Moreover, as we 
saw in the previous section, Jacobian matrices and their generalisation, the linear response, generate a 
graph structure that can be interpreted in causal terms. 

Dynamical effects. As a corollary in the increase of positive feedback loops it is observed that Hebbian 
synaptic adaptation leads to a systematic reduction of the dynamics complexity (transition from chaos 

19 If e is an edge, denote by o(e) the origin of the edge and t(e) its end. Then a feedback loop (or circuit) is a sequence 
of edges ei, ...,e^ such that o(ei-|_i) = t(ej), Vi = l...k — 1, and t(efc) = o(ei). Such a circuit is positive (negative) if the 
product of its edge's weight is positive (negative). 

20 In short, Hirsch [93] showed that cooperative systems, characterised by the property DFij(X.) > 0, Vi,j, VX g M, 
where DF is the Jacobian matrix, are convergent. As a matter of fact, this result holds for Jacobian matrices instead of 
synaptic weights matrix. But, in the particular example (132 1 ) the entry DFij of the Jacobian matrix is given by WijS' (Vj). 
Thus, it is proportional to Wij with a positive factor S'(Vj). 
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to fixed point by an inverse quasi-periodicity route, see fig. [5£i). As a corollary the largest Lyapunov 

(t) 

exponent X\ , which depends on the synaptic adaptation epoch r, decreases from positive to negative 
values. Two main effects contribute to this decay. The first effect is due to passive LTD term in (p0|) . 
The second one is related to an increase in the level of neurons' saturation. Basically cooperativity 
between neurons has a tendency to either render them more silent, or more active. In both case, this 
"pushes" them to the saturated part of the sigmoidal transfer function, reducing the average value of 
S'(Vi). Since the entry ij of the Jacobian matrix, DFij(V) = WijS'(Vi) an increase in the saturation 
of neurons has the effect of decreasing the spectral radius of DF(V) with a computable impact on the 
maximal Lyapunov exponent. 

Functional effects. This property has been exploited for pattern retrieval. Label by V the neuron 
state when the (time constant) input (external current) £ is applied to the network (see eq. (|32[) ) and by 
V the neuron state without The removal of £ modifies the attractor structure and the average value 
of observables. More precisely, let <j> be some suitable function and call: 

A W [^WV')) (TL WV)) W (61) 

where (0(V')) (r) is the (time or SRB) average value of 4> without £ and (</>(V)) M the average value in 
the presence of £. Two cases can arise. 

In the first case, the system is away from a bifurcation point and removal results in a variation of 
AM [4>] that remains proportional to provided £ is sufficiently small. In the present context, the linear 
response theory (see section [4| predicts that the variation of the average value of V is given b\l 21 l [43l [44] : 

AM [V] = - X M£ (62) 

where 

oc 

X M = £(£>F">M (63) 

71 = 

is a matrix whose entries are given by eq. ([57]) in section 0] Note therefore that AM [V] = — £ — AfM£ 
where the matrix A/M = Y^—i (-DF")M integrates dynamical effects. The application of £ implies a 
reorganisation of the dynamics which results in a complex formula for the variation of (V) M ; even if the 
dominant term is as expected. More precisely, as emphasised several times above, one remarks that 
each path in the sum , n \ is weighted by the product of a topological contribution depending only on 
the weights Wij and on a dynamical contribution. The weight of a path 7^ depends on the average value 

°f (riiLi f'( u ki-i {I ~ 1)))^ thus on correlations between the state of saturation of the units ko, . . . , k„-i 
at times 0, . . . , n — 1. 

In the second case the system is close to a bifurcation point and the presentation/removal of the 
input induces a sharp variation (bifurcation) in dynamics. In this case, eq. ((62]) does not apply anymore. 

21 We consider here the case (f>(V) = V for simplicity. For a general (diffcrcntiable) function tf>, the corresponding formula 
is IT5T1 : 

+ OO 

AM [<t>] = - (^F n V0> (T) £ 
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We expect therefore input presentation/removal to have a maximal effect close to bifurcations points. 
This can be revealed by studying the quantity A( T ) [<f>] in eq. (j6Tj) for the observable (j)(V) = S'(V), 
which measures the level of saturation of neuron i in the state Vi- Indeed, this quantity is maximal 
when Vi is in the central part of the sigmoidal transfer function, where neuron i is the most sensitive 
to small variations. Hence this quantity, called A' T '[S"], measures how neuron excitability is modified 
when the input is removed. The evolution of A^ T ^[S"] during learning following rule eq. (|60p is shown 
on fig. [8p (full lines) for two values of the passive forgetting rate A. A( T )[S"] is found to increase to 
a maximum at early learning epochs, while it vanishes afterwards. Interestingly, comparison with the 
decay of the leading eigenvalue /Ui (dotted lines) of the average Jacobian matrix (DF)^ shows that the 
maximal values of A( T )[S"] are obtained when is close to 1. This also corresponds to a vanishing of 
the maximal Lyapunov exponent. Hence, these numerical simulations confirm that sensitivity to input 
removal is maximal when the leading eigenvalue is close to 1. Therefore, "Hebb-like" learning drives 
the global dynamics in a region of the parameters space were sensitivity to the stimulus, manifested by 
a drastic change in the average of functions hence by a drastic change in the underlying dynamics, is 
maximal. This property may be crucial regarding memory properties of recurrent neural networks, which 
must be able to detect, through their collective response, whether a learnt input is present or absent. 
This property is obtained at the frontier where the strange attractor begins to destabilise (\fii\ = 1), 
hence at the so-called "edge of chaos" . 

Note that continuing adaptation after this phase of highest sensitivity, leads to a decay of sensitivity 
and to a stabilisation of the system to a fixed point (fig. [5^). This is not surprising. Insisting too much 
on adaptation to this stimulus drive the system to a state where activity pattern is essentially identical 
to the input, with no room any more for spontaneous activity. 

Conclusion. In fig. [Sfc we have represented the effect of presenting/removing the input before adap- 
tation, and after adaptation, in the region where the reactivity is maximal. Clearly, the presentation 
of £ after adaptation induces a sharp transition in dynamics, which is not occurring before adaptation. 
Moreover, this effect, inherited via learning, is robust to a small amount of noise, and selective (it does 
not occur for drastically different patterns) [59] . Though established in the context of a rather simple 
"neural" model, these results raise interesting questions and comments. They suggest that adaptation 
corresponds to some path in the space of parameters of (|32[) leading the system in a region were it acquires 
sensitivity to the input it has adapted to, this sensitivity being manifested by sharp and rapid variations 
in neurons dynamics. An obvious question is does this effect generalise to more complex (and realistic) 
architecture ? This is an ongoing research field. 

5.3 Effects of synaptic plasticity on spike trains statistics. 

Synaptic plasticity acts also on the statistics of spike trains. Let us now briefly mention recent works 
where the effects of plasticity on spike train statistics is analysed. 

Synapses update as an integration over spikes trains. Let us reconsider the equation ((9J for 
synaptic weights dynamics. The synaptic variation SWij is the integrated response of the synapse from 
neuron j to neuron i when neuron j sends a spike sequence [k>j] t _ T t and neuron i fires according to 
[ W i]+_T f This response is not a deterministic function, while @ is deterministic. As a matter of fact, 
the explicit form of g is usually derived from phenomcnological considerations as well as experimental 
results where synaptic changes can be induced by specific simulations conditions, defined through the 
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Figure 8: Fig. [8^. (left) Inverse quasi periodicity route induced by learning. The plotted quantity is m,( T )(to) = 
i S(Vi(to)), where to £ [0,T] is fixed and where the adaptation epoch index r varies. Fig. [Hp (middle) Sensitivity 

to the learned input. fj,\ is the largest eigenvalue of (DF)' 7 "' and L\ is the largest Lyapunov exponent. Fig. [St (right) 
Dynamical effect of input presentation before and after adaptation. 
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firing frequency of pre- and post-synaptic neurons [23J IM] > the membrane potential of the post-synaptic 
neuron [ID) , spike timing [1191 11241 [19] (see [123] for a review). Thus, these results are usually based on 
a repetition of experiments involving the excitation of pre- and post-synaptic neurons by specific spike 
trains. The phenomenological plasticity rules derived from these experiments are therefore of statistical 
nature. Namely, they do not tell us what will be the exact changes induced on synapses when this or this 
spike train is applied to pre- and post-synaptic neuron. Instead, they provide us the average synaptic 
change. Thus, the function g(Wij, [tOi] t _ T t , [u>j] t _ T t ) in fl9j is typically a statistical average of the 
synaptic response when the spike train of neuron j (resp. i) is [ujj] t _ T t (resp. [w»] t _ T t ), and the actual 
synaptic weight value is Wij. 

Slow synaptic update. On this basis let us assume that g is obtained via a time average 7r£ r - > (see 
eq. (|4"7j) in section [3]) of some function </> having a form depending on the type of "rule" considered (e.g. 
Hebbian or STDP). Namely, g has the form: 

9 (Wij, [w<] t _ T „ ( , [ui] t -T..t) = e ^ T) [<fa(Wj, ■)] ■ (64) 



where e is a parameter that will be typically small. In general 4nj can be expanded in terms of singlets, 
pairs, triplets, etc of spikes [771 HP] . 

Statistical effects of synaptic plasticity The coupled dynamics (]59[) has an impact on the set of 
admissible raster plots and the spikes train statistics. Typically, the empirical average constructed via 

( \ - (T) (T) 

the raster plot ui (T ' in the adaptation epoch r, changes from n ^ — > 7r Thus, the adaptation 

dynamics results in a sequence of empirical measures | an d the corresponding statistical model 

also evolves. Let us characterise this evolution using the tools introduced in section [3J The main idea is 

(T) 

to make the assumption that each 7r^ (T ) can be approximated by a Gibbs measure v^t) with potential 
ip( T \ Then synaptic adaptation writes: 



5WV = e V 



<kj(wjr\.) ■ (65) 



The synaptic update results in a change of parameters 7, ^ T+1S > = ^ T > +5^^ where S-y^ is assumed 
to be small (this is the role of the constant e in ([9])). This induces variation in statistical properties of 
raster plots (e.g. the topological pressure (|4"5)) -see [ID] for details). These variations can be smooth or 
not. 



Smooth variations If they are smooth one can show that there exist a function (W) such that the 
adaptation rule (|65p can be written in the form: 

SW (T) - eV w=w(T )4 T) (W). (66) 

Moreover, this quantity decay under synaptic adaptation. Thus, the adaptation rule (J66j is a gradient 
system where the function decreases when iterating synaptic adaptation rules. Were the transition 
t — > t + 1 to be smooth for all r, would reach a minimu n@ at some W* as r — > 00. Such a minimum 
corresponds to Vw-T^ = 0, thus to 

22 In implementing synaptic update rule, one adds conditions ensuring that weights do not diverge. This condition ensures 
that is bounded from below, W 
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SW*j = u+* [fa] = 0, ,V*,i = 1...JV, 
according to eq. (|66p . Hence, this minimum corresponds to a static distribution for the synaptic weights. 

Static synaptic weights. Since this imposes a condition on the average value of the 4>ij 's, this imposes 
as well the statistical model as a Gibbs measure v^, with a potential: 

t/>* = $ + A*.0, (67) 

where t/>* = (V J ij)^-_ 1 - The potential $ in (f6"T)) is such that = is u> is admissible and = — oo 

if it is forbidden, so that forbidden raster plots have zero probability. We use the notation = (0y )f^ =1 , 

A* = (Ay) . ._ and A*.</> = X)fj=i Kj't'ij- The statistical parameters Ay, are given by eq. (|50[) in section 
1331 

As a conclusion, the statistical model, in the sense of section RT2I is a Gibbs distribution such that 
the probability of a spin block R of depth n obeys eq. (|52l) with a potential (see eq. (|7T|) for an 
explicit form.)) When this corresponds to the asymptotic state for a synaptic adaptation process, this 
potential provides us the form of the statistical model after adaptation, and integrates all past changes 
in the synaptic weights. 

Moreover, it has a deep implication. Since the Gibbs distribution obeys the variational principle (|48[) 
with v^,' [4>] — , the probability distribution has maximal entropy. In other words, synaptic adap- 
tation rules of type (|65p . when they converge, drive the system toward a dynamics where the statistical 
entropy of spike train is maximal, taking into account the constraints^] imposed by dynamics. 

Singular variations. The synaptic weights variations can induce a change in the set of admissible 
raster plots that the system is able to display. When this happens the set of admissible raster plots is 
suddenly modified by the synaptic adaptation mechanism. Formerly forbidden sequences become allowed, 
formerly allowed sequences become forbidden, but also a large core of legal sequences may remain legal. 
These changes depend obviously on the detailed form of neuron dynamics (|14[) and of the synaptic update 
mechanism An interesting situation occurs when the set of admissible raster plots obtained after 
adaptation belongs to fl S^c^+i). In this case, adaptation plays the role of a selective mechanism 

where the set of admissible raster plots, viewed as a neural code, is gradually reducing, producing after n 
steps of adaptation a set n" i=1 S 7 ( m ) which can be rather small. If we consider the situation where (fT4"]) 
is a neural network submitted to some stimulus, where a raster plot uj encodes the spike response to the 
stimulus, then Y,^ is the set of all possible raster plots encoding this stimulus. Adaptation results in a 
reduction of the possible coding, thus reducing the variability in the possible responses. 

Example: Spike Time Dependent Plasticity. As an example we consider a neural network of 
type (|4"Tj) with an adaptation rule inspired from (TTTj) with an additional term rdW^j\ —1 < r^ < 0, 
corresponding to passive LTD. 

23 We mean that, as emphasised several times in the paper, dynamics is not able to produce all possible spikes trains. 
Henceforth, statistical entropy must be maximised on a subset of spike trains which can be relatively small compared to 
the whole set of possible spike trains (2 NT possible spike trains of length T for a system of TV neurons). Note also that 
"maximal entropy" does not mean "equi-probability" here. 
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SW^=e 



r d W 1 



(r) 



T+T 



^E^ T) (*) E /("K (T) (* 



t=T. 



(68) 



where /(x) is given by (jX2j) and with: 



Set 



clef 



T s = 2 max(r + , t_) . 



then {5HK has the form Qg), tfwi? = ettJJ) [&i(Wtf> ■)]■ 



(69) 



u=-T s 



Static weights. Thanks to the soft bound term r^Wij the synaptic adaptation rule admits a static 
solution given by: 



Wij = - 



TI=-t s /("Vjw [uj(0)ui(u 
rd 



(70) 



The sign of Wij depend on the parameters A+, T s , but also on the relative strength of the terms 



(T) 

77 m Note that this equation may have several solutions 



Spike train statistics in a static weights regime. As emphasised in section 13.21 and 15. 3[ when 
the synaptic adaptation rule converges to a fixed point, the corresponding statistical model is a Gibbs 
measure with a potential 

N N 

= EE A *^> 

i=l 3 = 1 

where the value AL of the Lagrange multipliers is constrained by the relation ([70)1 . Henceforth, the 
probability of an admissible spike bloc R is given by: 



P\R\S] 



exp 



N N Ti-l T s 

EEWE E f{u)^{t)^{u+t) 

i=l j=l t=0 u=-T s 



(71) 



Note that the term r d Wij arising in the definition of the potential can be removed thanks to the normal- 
isation constant Z n [X±(S), . . . , A; (5)]. 



Conclusion At the end of section [3] we were asking whether it is possible to propose a canonical form 
for spike trains statistical models relying on some generic principle ? Here, we have exhibited an example 
where such a construction can be made. Our argumentation suggests that Gibbs measures may be good 
statistical models for a neuron dynamics resulting from slow adaptation rules where synaptic weights 
converge to a fix value. In this case, synaptic weights contains the whole history of the neural network, 
expressed in the structure of the generating function of cumulants (the topological pressure) and in the 
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structure of allowed/forbidden raster plots (potential $ in (|57|) ). This theoretical results, confirmed by 
numerical experiments [40] . can be compared with the recent paper of Schncidman and collaborators, 
already quoted in section [31 proposing a Gibbs distribution with Ising like potential to match empirical 
data on the salamander retina |155j . Our approach suggests that Gibbs distribution could be ubiquitous 
and that the corresponding potential can be inferred according to the plasticity mechanisms at work, 
defining the "rule" . 

6 Conclusion 

In this paper we have presented a series of results related to questions, that, we believe, are central in the 
computational neuroscience community. These questions, when addressed from the point of dynamical 
systems theory, shed new light on neural network dynamics with possible outcomes towards experimen- 
tation. Interestingly enough, the tools and concepts used here however are not restricted to the field of 
neural networks but could be applied to other type of "complex" systems. 

We now would like now to point out some "challenging" points raised in this paper, which, to our 
opinion, are central at the current state of the art, in the neuroscience community, especially for those 
people who want to use models and their analysis to "understand some fundamental keys ruling the be- 
haviour of neural networks" . We also believe that some of these questions address also to experimentalists 
community 

Finite size corrections. When dealing with large populations of networks, theoretical methods such 
as mean-field approaches, neural mass models, large deviations, use a limit N — > oo where N is the 
number of neurons. When dealing with "concrete" neural systems the number is finite and finite size 
systems can have a behaviour that departs strongly from the limiting system. Computing these finite 
size corrections is an open problem in the whole scientific community (not only neuroscience). 

Finite time corrections. A similar question holds for finite time effects. While many theoretical 
methods assume stationarity in the data, concrete experiments handle non stationary or transient dy- 
namics. There exists currently empirical methods to tackle this problem such as sliding time windows. 
But, to the best of our knowledge there are rather few methods (i) to estimate the width of this window 
which must be large enough to ensure reliable statistical estimates and smaller than relevant character- 
istic time scales of the non-stationary dynamics; (ii) to define a priori the statistical models (resp. the 
set of observable) used to characterize dynamics; (iii) to propose and compute reliable indicators that 
guarantee the liability of the result. 

Correlated weights. Synaptic weights are non independent in real neuronal networks since e.g. plas- 
ticity build correlations, that can have long range in space and time. How to characterize these cor- 
relations ? How to measure them ? How to handle them in a model ? For example, the mean-field 
approaches developed in section cannot be extended in a straightforward way to this case. Here again 
there is a need to invent new theoretical methods. Promising results using spectral properties of graphs 
[1081 IT21 11071 [22] combined to such linear responses methods as developed in section 0] could provide 
a breakthrough. 
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Continuous time. Most phcnomcnological models of neurons, such as Hodgkin- Huxley's, use a descrip- 
tion with ordinary differential equations assuming a continuous time. This comes from the description 
of neurons in terms of membrane potential, a notion which corresponds to an integration over space dx 
and time scale dt where classical equations of electrodynamics holds (Kirchhoff and Ohm's law). Also, as 
we saw, Hodgkin-Huxley's equations for example, corresponds to integrating the activity of ionic channel 
at space and time scales so that the notion of probability of open/closed channel have a meaning and, 
moreover, so that the Markovian approach used for the master equations (|2I3I40 in section [TT2l holds. In 
other words, writing neurons dynamics at the scale of ionic channel would requires the use of a different 
physics. Though this remarks looks obvious it prevents one to take the continuum limits dx — > or 
dt — ► without caution. On the opposite, the "spike" description of neurons is discrete and assume a 
time discretisation S (see section [l~2l such that a given neuron can produce, at most, one spike within this 
time delay. This time discretisation is essential for the definition of raster plots, as we saw. The existence 
of a minimal time scale S can be defended using arguments from physics and neuronal characterisation 
(see [45]). However, time discretisation as done in section ll~2l for the definition of raster plots, imposes a 
time grid which can be contested as well since one may argue that a neuron spike can occur at "any" time 
(with the restriction discussed above on "continuous" time) while the time grid set it to a time multiple 
of 5 [114] . So, a central question is: does the effects induced by this slight error on a spike time matters ? 
This question can be addressed in the realm of models (like gIF models) and the answer is: "it depends" . 
Indeed, according to the control parameters values such an error can be damped (then it is harmless) or 
amplified by dynamics. This is precisely the discussion on contraction and initial conditions sensitivity 
presented in section [2.2.31 Here, the mathematical analysis relies strongly on the simple structure of IF 
models. But the question of time discretisation and "How precise is the timing of action potentials" |114j 
must be certainly addressed in a more general context and for more general models. 

Neural code and predictability Another interesting issue, raised by the dynamical system point of 
view is the following. When building models which reproduce the dynamics of neural assemblies, one 
is faced to the question "how well does this model approximate real neural systems" . Let us state in a 
different way. Assume that we have built some artificial neural network that mimics some part of the 
brain and assume (Gedanken experiment) that we are able to remove this part of the brain and to replace 
it by our artificial system, what do we need to ensure that this "cyborg brain" works "as" the original 
one ? This is a (too) wide open question but let us focus the discussion on "spikes" aspects for simplicity. 
Our artificial system receives spikes from other parts of the brain and respond with spikes trains that it 
sends to various brain areas. Must this device be able to reproduce exactly, spikes par spikes, the response 
of the piece of brain that it replaces ? It is possible to reproduce exactly the response of one neuron to 
Poisson stimuli [106] . and also at the level of a network, to reproduce exactly finite spike trains coming 
from biological neurons over a finite time |147j . But, from the dynamical system point of view this 
appears to be too restrictive since, as we saw, in many cases dynamics is chaotic. If one jitters the time of 
a single spike, the output of the network can change dramatically. Thus this property is not robust and 
reliable. So which characteristics of the spike train must be reproduced ? This question seems a clear 
challenge for a not so near future. This also raise questions such as: Are all details important in mod- 
elling neural networks ? Which details can be neglected and what is imposed by the biological structure ? 

Though the analysis presented here is rather simple compared to the overwhelming richness of brain 
dynamics, we hope that this work will be useful for readers interested in this beautiful ongoing research 
field, computational neuroscience. 
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